Spin Change of Asteroid 2012 TC4 Probably by Radiation Torques

In the first step, we analyzed the photometry data from 2012 and 2017 using the two-period Fourier series method (Pravec et al. 2005, 2014). Concerning the 2012 observations, we used all but one photometric light-curve series taken from 2012 October 9.9 to 2012 October 11.2 (see Table 1). In particular, we excluded the data taken on 2012 October 10.7 by Polishook (2013), where we found a possible timing problem.¹² Concerning the 2017 observations, we used a selected subset of the data taken between 2017 October 9.1 and 2017 October 11.1. This choice was motivated by noting that the observing geometry during this time interval in 2017 was very similar (due to the fortuitous resonant return of the asteroid to Earth after 5 yr) to the geometry of the time interval of the 2012 observations. This choice minimizes the possible (but anyway small) systematic effects due to changes in observing geometry when comparing the results of our analysis of the data from the two apparitions (see below).

We reduced the data to the unit geo- and heliocentric distances and a consistent solar phase using the H–G phase relation assuming G = 0.24, converted them to flux units, and fitted them with the fourth-order two-period Fourier series. A search for periods quickly converged, and we found the two main periods P₁ = 12.2183 ± 0.0002 and P₂ = 8.4944 ± 0.0002 minutes in 2012 and P₁ = 12.2492 ± 0.0001 and P₂ = 8.4752 ± 0.0001 minutes in 2017. The phased data and best-fit Fourier series, together with the postfit residuals, are plotted in Figure 1. We note that the smaller formal errors of the periods determined from the 2017 data were due to the higher quality of the 2017 observations (the best-fit rms residuals were 0.081 and 0.065 mag for the 2012 and 2017 data, respectively; see also the postfit residuals plotted in the bottom part of Figure 1). As for possible systematic errors of the determined periods, the largest could be due to the so-called synodic effect. It is caused by the change of position of a studied asteroid with respect to the Earth and Sun in the inertial frame during the observational time interval. An estimate of the magnitude of the synodic effect can be obtained using the phase-angle-bisector approximation, for which we used Equation (4) from Pravec et al. (1996). Using this approach, we estimate the systematic errors of the determined periods ΔP₁ = 0.0005 and 0.0002 minutes and ΔP₂ = 0.0002 and 0.0001 minutes in 2012 and 2017, respectively. These systematic errors are only slightly larger than the formal errors given above. A caveat is that the formula of Equation (4) in Pravec et al. (1996) was determined for the case of a principal-axis rotator. An exact estimate of the systematic period uncertainties due to the synodic effect for a tumbling asteroid would require an analysis of its actual non-principal-axis rotation in a given observing geometry, but we did not pursue it here, as the effect was naturally surmounted by the physical modeling presented in the next section.

Zoom In Zoom Out Reset image size

As will be shown in the next section, the strongest observed frequency in the light curve, ${P}_{1}^{-1}$, is actually a difference between the precession and the rotation frequency of the tumbler: ${P}_{1}^{-1}={P}_{\phi }^{-1}-{P}_{\psi }^{-1}$. The second-strongest frequency, ${P}_{2}^{-1}$, is then the precession frequency, ${P}_{\phi }^{-1}$. This is a characteristic feature of the light curve of a tumbling asteroid in short-axis mode (SAM; see below). We note that the same behavior was observed for (99942) Apophis and (5247) Krylov, which are also in SAM (Pravec et al. 2014; Lee et al. 2020).

The principal light-curve periods P₁ and P₂ determined from the 2012 and 2017 observations differ at a high level of significance, formally on the order of about 100σ. While the systematic errors due to the synodic effect could decrease the formal significance by a factor of a few, the significance of the observed period changes would still remain large, on the order of several tens of σ. To interpret these findings in more depth, we turned to construct a physical model of 2012 TC4 as a tumbling object in the next section.

The light-curve data set from 2017 is much richer than that from 2012, so we started with the inversion of the 2017 data. We investigated possible frequency combinations based on the fact that the main frequencies f₁ and f₂ of a tumbling asteroid light curve are usually found at 2f_ϕ and 2(f_ϕ ± f_ψ) or low harmonics and combination, where f_ϕ is the precession and f_ψ is the rotation frequency, respectively, and the plus sign is for long-axis mode (LAM) and the minus sign for SAM (Kaasalainen 2001). Using the values f₁ = 4.898 and f₂ = 7.079 hr⁻¹ from the previous section, we found eight possible frequency combinations: f₁ = f_ϕ, f₂ = 2(f_ϕ − f_ψ) (SAM1); f₁ = 2(f_ϕ − f_ψ), f₂ = 2f_ϕ (SAM2); f₁ = 2(f_ϕ − f_ψ), f₂ = f_ϕ (SAM3); f₁ = f_ϕ − f_ψ, f₂ = f_ϕ (SAM4); f₁ = 2f_ϕ, f₂ = 2(f_ϕ + f_ψ) (LAM1); f₁ = 2f_ϕ, f₂ = f_ϕ + f_ψ (LAM2); f₁ = f_ϕ, f₂ = (f_ϕ + f_ψ) (LAM3); and f₁ = f_ϕ + f_ψ, f₂ = f_ϕ (LAM4). Then we conducted the shape and spin optimization for these combinations in the same way as in Lee et al. (2020). It was found that only the SAM4 solution provided an acceptable fit to the data and was physically self-consistent.

We used 34 light curves from 2017 October. Four light curves from 2017 September were very noisy and did not further constrain the model, so we did not include them in our analysis. We inverted the light curves with the method and code developed by Kaasalainen (2001) combined with Hapke's light-scattering model (Hapke 1993). According to Reddy et al. (2019), the colors of TC4 are consistent with the C or X complex; also, its spectrum is X type, with Xc type being the best match. The X complex contains low- and high-albedo objects, but the E type seems most likely because the high circular polarization ratio suggests that TC4 is optically bright (Reddy et al. 2019). This is in agreement with Urakawa et al. (2019), who also reported X-type colors. As a result, we used Hapke's model with parameters derived for E-type asteroid (2867) Šteins (Spjuth et al. 2012): ϖ = 0.57, g = − 0.30, h = 0.062, B₀ = 0.6, and $\bar{\theta }=28^\circ$. Because our data did not cover low solar phase angles, the h and B₀ parameters for the opposition surge and also roughness $\bar{\theta }$ were fixed. We optimized the ϖ and g parameters, and they converged to values of ϖ = 0.69, g = − 0.20, which gave a geometric albedo of 0.34. An alternative solution with fixed values ϖ = 0.57, g = − 0.30 provided an only marginally worse fit, and the kinematic parameters were not affected. In general, the solution of our inverse problem was not sensitive to Hapke's parameters like in previous studies (e.g., Scheirich et al. 2010; Pravec et al. 2014; Lee et al. 2020).

The rotation and precession periods had the following values: P_ψ = 27.5070 ± 0.002 and P_ϕ = 8.47512 ± 0.0002 minutes, respectively. The 1σ uncertainties were estimated from the increase of χ² when varying the solved-for parameters; given a number of measurements of about 8600, the 3σ uncertainty interval corresponds to an about 5% increase in χ² (e.g., Vokrouhlický et al. 2017). The direction of the angular momentum vector in ecliptic coordinates was λ = 92°, β = − 89.6°, practically oriented toward the south ecliptic pole. Normalized moments of inertia were I₁ = 0.41, I₂ = 0.81, but the model was not too sensitive to their particular values. The inertia moments computed from the 3D shape (assuming constant density) were I₁ = 0.435, I₂ = 0.831, indicating consistency with the kinematic parameters above.

The dark facet, which is always introduced into a convex shape model to regularize the solution (Kaasalainen & Torppa 2001), represented about a few percent of the total surface area, and forcing it to smaller values led to a worse fit. This might mean that there is some albedo variation on the surface of TC4 or that its real shape is highly nonconvex and a convex-shape approximation has its limits.

For this model, we used 14 light curves from 2012. The rotation and precession periods were now P_ψ = 27.873 ± 0.005 and P_ϕ = 8.4945 ± 0.0003 minutes, respectively. The 1σ uncertainties were estimated in the same way as for the 2017 model, namely, from the increase of χ². The direction of the angular momentum vector was λ = 89°, β = − 90°, with moments of inertia I₁ = 0.43, I₂ = 0.80. The 3D shape model was similar to that reconstructed from 2017 data, and the direction of the angular momentum vector was practically the same.

In spite of consistency in all other solved-for parameters, we thus observed that the periods P_ϕ and P_ψ for the two apparitions were significantly different. Attempts to use the 2017 values to fit the 2012 data, or vice versa, led to a dramatically worse fit. We scanned the period parameter space around the best-fit values and plotted the relative χ² values; see Figure 2. The minima in χ² for the 2012 and 2017 periods are clearly separated.

Zoom In Zoom Out Reset image size

The 3D models reconstructed from the two apparitions independently have similar global shapes, but their details are different (Figure 3). If we take the shape reconstructed from the 2017 data and use it to fit the 2012 data, it gives a satisfactory fit to the light curves. However, a better way to use all data together is not to treat them separately but to instead invert both apparitions together with a common shape model that would differ only in kinematic parameters. Therefore, we modified the original inversion code of Kaasalainen (2001) to enable including two independent light-curve sets. We assumed that the only parameters that were different for the two apparitions were the rotation and precession periods P_ψ, P_ϕ and initial Euler angles ϕ₀, ψ₀. The parameters describing the shape and the direction of the angular momentum vector were the same. So the full set of kinematic parameters for a two-apparition model was (λ, β, ${\phi }_{0}^{(1)}$, ${\phi }_{0}^{(2)}$, ${\psi }_{0}^{(1)}$, ${\psi }_{0}^{(2)}$, ${P}_{\psi }^{(1)}$, ${P}_{\psi }^{(2)}$, ${P}_{\phi }^{(1)}$, ${P}_{\phi }^{(2)},{I}_{1},{I}_{2})$, where the superscript (1) is for the 2012 apparition and (2) is for the 2017 apparition. The two light-curve data sets were independent in the sense that the integration of kinematic equations (Equation (A.3) in Kaasalainen 2001) was done separately for 2012 and 2017 data; the two epochs were not directly connected. The final model shape is almost identical to that reconstructed from only 2017 data shown in Figure 3. The fit of the final model to individual light curves is shown in Appendix B. Rotation and precession periods converged to practically the same values as with the independent treatment of each of the two apparitions. The best-fit parameters for the 2012 and 2017 apparitions are listed in Table 2. The physical models of 2012 TC4 from 2012 and 2017 and the light-curve data set used to reconstruct these models are available from the DAMIT database¹³ (Ďurech et al. 2010).

Zoom In Zoom Out Reset image size

Table 2.  Parameters of the Model in Figure 3 Reconstructed from 2012 and 2017 Light Curves

	2012	2017
Pψ (minutes)	27.8720 ± 0.0007	27.5070 ± 0.0002
Pϕ (minutes)	8.4944 ± 0.0005	8.47511 ± 0.00008
ϕ0 (deg)	322 ± 8	74 ± 9
ψ0 (deg)	198 ± 0.02	180 ± 0.2
JD0	2,456,210.00	2,458,032.69
λ (deg)	103 ± 78
β (deg)	−88.5 ± 0.7
I1	0.42
I2	0.81
w	0.67
g	−0.20
h	0.062
B0	0.6
$\bar{\theta }$ (deg)	28
δL/L	0.00078 ± 0.00006
δE/E	−0.0035 ± 0.0001

Note. The reported errors correspond to standard deviations of parameters computed from bootstrap models. Parameters without errors were fixed (Hapke's parameters) or did not change (I₁ and I₂). The formally small uncertainty of ψ₀ means that this initial orientation angle is correlated with the shape and does not change significantly with different bootstrap data sets. Parameter λ is practically unconstrained because the angular momentum direction is very close to β = − 90°. Values δL/L and δE/E are relative changes of angular momentum L and energy E between 2012 and 2017, i.e., δL/L = (L₂₀₁₇ − L₂₀₁₂)/L₂₀₁₇ and δE/E = (E₂₀₁₇ − E₂₀₁₂)/E₂₀₁₇.

Download table as:  ASCIITypeset image

To estimate the uncertainties of physical periods and further robustly demonstrate that their change between 2012 and 2017 is significant, we created bootstrapped data samples and repeated the light-curve inversion for them. From both the 2012 and 2017 data sets, we created 1000 bootstrap samples by randomly selecting the same number of light curves from the original data set. For the 2017 October bootstrap, 279 final shape models had a clearly wrong inertia tensor that was not consistent with the kinematic I₁, I₂ parameters, so we removed them from the analysis. For all remaining bootstrap models, we plot histograms of P_ϕ and P_ψ distribution in Figure 4. The standard deviations of the precession period P_ϕ are 0.0005 and 0.00009 minutes for the 2012 and 2017 data, respectively, which are similar to the uncertainty values derived in Sections 3.1 and 3.2. For the rotation period P_ψ, these standard deviations are 0.0008 and 0.0002 minutes, which is significantly smaller than our χ²-based estimate. Nevertheless, the difference between periods determined from the 2012 and 2017 apparitions is much larger than their uncertainty intervals, and we are not aware of any random or model errors that could cause such a difference. Our conclusion is that the spin state of TC4 has changed from 2012 to 2017, and the rotation and precession periods have decreased. In what follows, we try to interpret this change using a theoretical model of TC4 spin evolution with the relevant torques.

Zoom In Zoom Out Reset image size

The unique observational opportunities of 2012 TC4 are directly related to its exceptional orbit. The asteroid had a deep encounter with Earth on 2012 October 12, during which the closest distance to the geocenter was approximately 95,000 km (Figure 5). However, the more unusual circumstance was that the 2012 close encounter resulted in a change of 2012 TC4's orbital semimajor axis that placed it nearly exactly at the 5:3 resonance with Earth's heliocentric motion. As a result, 5 yr after the first close encounter, i.e., on 2017 October 12, the relative configuration of the asteroid and Earth nearly exactly repeated, again placing it in a deep encounter configuration. This time, the closest approach to the geocenter had an even closer distance of 50,200 km. Astrometric observations during the two close approaches, including Arecibo and Green Bank radar data taken in 2017, allowed a very accurate orbital solution over the 5 yr period of time in between 2012 and 2017. In the context of this paper, we note that it also provided interesting information about the nongravitational effects that need to be empirically included in orbit determination. Adopting a methodology from cometary motion (e.g., Marsden et al. 1973 and Farnocchia et al. 2013 or Mommert et al. 2014 for the asteroidal context), we note the following values of radial A₁ and transverse A₂ accelerations: (i) A₁ = (2.17 ± 0.80) × 10⁻¹¹ au day⁻², or (4.35 ± 1.60) × 10⁻¹⁰ m s⁻², and (ii) A₂ = − (2.73 ± 0.65) × 10⁻¹³ au day⁻² (both assume ∝ r⁻² heliocentric decrease; see the JPL/Horizons webpage, https://ssd.jpl.nasa.gov/sbdb.cgi). At first sight, these values appear very reasonable (compare with, e.g., A₁ and A₂ fits for the ≃4 m body 2009 BD; Mommert et al. 2014; Vokrouhlický et al. 2015a).

Zoom In Zoom Out Reset image size

If we were to interpret both components as a result of radiation forces, we might further obtain useful information about the body. The radial component would represent the direct solar radiation pressure. In a simple model, where we assume a spherical body of size D and bulk density ρ, we have A₁ ≃ 3C_RF₀/(2ρDc), with C_R the radiation pressure coefficient (dependent on sunlight scattering properties on the surface), F₀ ≃ 1367 W m⁻² the solar constant, and c the light velocity. Adopting C_R ≃ 1.2 and D ≃ 10 m, we obtain a very reasonable bulk density ρ ≃ (1.9 ± 0.7) g cm⁻³. The transverse component of nongravitational orbital effects makes sense when interpreted as the Yarkovsky effect (e.g., Vokrouhlický et al. 2015a). Note that the abovementioned value of A₂ translates to a secular change of the semimajor axis da/dt = −(110 ± 26) × 10⁻⁴ au Myr⁻¹ (see, e.g., Farnocchia et al. 2013). Given the near extreme obliquity of the rotational angular momentum, we may safely restrict to the diurnal component of the Yarkovsky effect. The negative value of A₂, or da/dt, corresponds well to the retrograde sense of 2012 TC4 rotation (implied by the direction of the rotational angular momentum vector; Section 3). Next, borrowing the simple model for a spherical body from Vokrouhlický (1998) and fixing the size D = 10 m, the inferred bulk density ρ ≃ 1.9 g cm⁻³, and the surface thermal conductivity K = 0.05 W m⁻¹ K⁻¹, we estimate that the corresponding surface thermal inertia is ${\rm{\Gamma }}\simeq {490}_{-250}^{+270}$ in SI units (though we note that there is also a lower-inertia solution possible, like in Mommert et al. 2014). This is a very adequate value too (e.g., Delbó et al. 2015). Obviously, due to the many frozen parameters in the model (and its simplicity), the realistic uncertainty in Γ would be larger. However, it is not our intention to fully solve this problem. We satisfy ourselves with the observation that the needed empirical nongravitational accelerations in the orbital fit may be very satisfactorily interpreted as radiation-related effects. In the next sections, we show that the change in the rotation state between the 2012 and 2017 observation epochs may be very well explained by the radiation torque known as the YORP effect (e.g., Vokrouhlický et al. 2015a).

Our methods and mathematical approach used to describe the evolution of the rotation state of 2012 TC4 are presented in Appendix A; therefore, here we provide just a general outline. We numerically integrated Euler Equations (A2) and (A4) describing the spin state evolution of 2012 TC4 in between the observation runs in 2012 and 2017. The kinematical part, describing the transformation between the inertial frame and the body frame defined by the principal axes of the tensor of inertia, was parameterized by the Rodrigues–Hamilton parameters λ = (λ₀, λ₁, λ₂, λ₃) (see, e.g., Whittaker 1917). This choice helps to remove problems related to coordinate singularity given by the zero value of the nutation angle. The dynamical part is represented by the evolution of the angular velocity ω in the body frame. Note that in most cases of asteroid light-curve interpretation, a simplified model of a free top would be sufficient (also used in Section 3 to fit the 2012 and 2017 data separately). However, the evidence of the change of the rotation state of 2012 TC4 in between the 2012 and 2017 epochs, discussed above, requires appropriate torques to be included in the model. We addressed two effects:

• 1.  
  gravitational torques due to the Sun and Earth and
• 2.  
  radiation torques due to the sunlight scattered by the surface and thermally reradiated (the YORP effect; e.g., Bottke et al. 2006; Vokrouhlický et al. 2015a).

The tidal gravitational fields of the Sun and the Earth were represented in the body frame using the quadrupole approximation (Equation (A7); e.g., Fitzpatrick 1970; Takahashi et al. 2013). The nature of the perturbation is different for the Sun and Earth. In the solar case, the gravitational torque results in a small tilt of less than 1' describing a small segment on the precession cone. The effect of the Earth-induced gravitational torques manifests as an impulsive effect only during close encounters. Our main goal was to verify that, due to the fast rotation of 2012 TC4, the effect averages out and cannot contribute to the observed change in rotational frequencies. Indeed, Figure 6 shows the change of the osculating rotational (intrinsic) angular momentum and energy during the 2012 close encounter with Earth (only about six times larger effects are observed during the closer encounter in 2017, but this is not relevant for our analysis anyway, because the observations preceded this approach). The initial data of the simulation were taken from the observations fit in 2012, namely, before the encounter. Recall that the characteristic timescale of the encounter is about half a day (determined, for instance, as a width of the characteristic Earth relative velocity increase in the bottom panel of Figure 5 at an ≃6.8 km s⁻¹ level, half the value between the asymptotic and peak velocities). In comparison, the characteristic rotational periods P_ψ and P_ϕ are of the order of minutes, i.e., much shorter. As a result, the effect of Earth's gravitational torque efficiently averages out during each of the rotation cycles (a significant effect may be expected only for very slowly rotating bodies and sufficiently deep encounters, such as seen for 4179 Toutatis during its 2004 close encounter with Earth; e.g., Takahashi et al. 2013). Therefore, we may conclude that the gravitational torques cannot explain the observed change in the rotation state of 2012 TC4 in between the 2012 and 2017 epochs. Still, we keep them in our model for the sake of completeness.

Zoom In Zoom Out Reset image size

The radiation torques are of a quite different importance. It is well known that the YORP effect is able to secularly change the rotational frequency and tilt the rotational angular momentum in space. While mostly studied in the limit of a rotation about the shortest axis of the inertial tensor, generalizations to the tumbling situation were also developed. Both numerical (Vokrouhlický et al. 2007) and analytical (Cicalò & Scheeres 2010; Breiter et al. 2011) studies confirmed that the YORP effect is able to change rotational angular momentum, and its orientation, in both the inertial space and bodyframe in an appreciable manner. As usual, the effect is more important on small bodies. Here we included the simple, zero-inertia limit developed in Rubincam (2000) and Vokrouhlický & Čapek (2002); see Equation (A8). A first look into the importance of the radiation torques may be obtained by taking the nominal solution of the rotation state from the 2017 observations, including the appropriate shape model, and propagating it backward in time to early 2012 October (the epoch of the first set of observations). We use the more accurate 2017 model as a reference, rather than the one constructed from the poorer 2012 observations. The length scale of the model was adjusted to correspond to an equivalent sphere of diameter D = 10 m, and the density was ρ = 1.4 g cm⁻³. We verified that the results have the expected invariance to rescaling of both ρ and D such that ρ D² = const. Our nominal choice of the 1.4 g cm⁻³ bulk density, albeit conflicting with the suggested E-type spectral classification of TC4, is therefore linked to the assumed equivalent size of 10 m, but it might be redefined according to the rescaling principle. This combination of parameters provides a very nice match of the P_ϕ period change due to our YORP model (see Section 4.3).

Figure 7 provides information about the secular change in the rotational angular momentum L (top) and energy E (bottom) in that simulation. Here we see a long-term change in both quantities. The wavy pattern is due to the eccentricity of the 2012 TC4 orbit and a stronger YORP torque at perihelion. The accumulated fractional change in both L and E is a few times 10⁻³. This is promising because the observed change in these quantities is of the same order of magnitude (see Table 2) and makes us believe that the change in the directly observable P_ψ and P_ϕ periods will also be as needed. Nevertheless, we also note a difference. The simulation results shown in Figure 7 indicate that both angular momentum and energy increased from 2012 to 2017. Such behavior is perhaps expected in the first place (for instance, in the case of a body in a principal-rotation state, YORP would necessarily affect both E and L in the same way). However, rotation state solutions from observations in 2012 and 2017 tell us something else (see Table 2): the rotational angular momentum L increased between 2012 and 2017, while the energy E decreased. Before commenting more on this difference, we first provide a more detailed analysis of the radiation torque effects for 2012 TC4 in our modeling, this time using the whole suite of acceptable initial data and shape models (all compatible with the observations; Section 3.4). This will allow a statistical assessment of the predicted values.

Zoom In Zoom Out Reset image size

An ideal procedure of proving that the observed changes in tumbling-state periods P_ϕ and P_ψ are due to the radiation (YORP) torques would require a highly reliable theoretical model (numerical propagation of 2012 TC4's rotation state with appropriate torques included) employed to fit all available observations (in our case, data from 2012 and 2017 October). Obviously, the only "comfort" of this analysis would be to possibly adjust some free (unknown) parameters. Unfortunately, such a plan is presently too ambitious; thus, we resort to a simpler way.

Recall the much easier situation when the YORP effect has been searched (and detected) for asteroids rotating in the lowest-energy mode, namely, about the shortest principal axis of the inertia tensor. In this case, YORP results in a secular change of the unique rotation period P (as in P_ϕ and P_ψ, when the body tumbles). The measurements are rarely precise enough, and the effect strong enough, to directly reveal the change in the period P (see, though, an exception for 54509 YORP; Lowry et al. 2007). More often, one uses the fact that the linear-in-time change in P produces a quadratic-in-time effect in the rotation phase. Properly linking the asteroid rotation phase over many observation sessions with an empirical quadratic term helps to characterize changes of P that are individually too small to be determined from one-apparition observations (see, e.g., Vokrouhlický et al. 2015a). This approach adopts an empirical magnitude of the quadratic term in the rotation phase and simply solves it as a free parameter (not combining it with a theoretical model at that stage). Interpretation in terms of the YORP effect is done only a posteriori, when the fitted amplitude of the phase-quadratic term is compared with a prediction from the YORP model. And even then, the comparison is often not simple, because the model prediction for YORP is known to depend on unresolved small-scale irregularities of the body shape. On several occasions, one has to be satisfied with a factor of a few difference accounting for the model inaccuracy.

Things are quite a bit more complex when the body tumbles. First of all, it is not clear how to set the empirical approach from above and apply it in this situation. At the same time, the direct modeling approach is probably even less accurate than in the case of rotation about the principal axis of the inertia tensor. Not only does the worry about the role of unresolved small-scale irregularities remain, but the present YORP model is restricted to the zero thermal conductivity limit (see, e.g., Vokrouhlický et al. 2007 for the numerical approach and Cicalò & Scheeres 2010 and Breiter et al. 2011 for analytical studies). The Yarkovsky acceleration for tumblers was evaluated with thermal models (e.g., Vokrouhlický et al. 2005, 2015b), but in these cases, P_ϕ and P_ψ were slightly tweaked to make them resonant (an approach we cannot afford here). In this situation, we adopted the following simple procedure.

4.3.1. Model Based on 2017 Data

We start with a set of models uniquely based on the most reliable and accurate observations from the 2017 apparition. In particular, we constructed 687 variants of the 2012 TC4 physical model (Section 3.1). They are all very similar because they sample tight parameter variations, all resulting in acceptable fits of the data. These represent (i) slightly modified initial rotation parameters (Euler angles and their derivatives, as well as the inertial space direction of the rotational angular momentum vector) and (ii) slight shape variants of the body. The initial epoch MJD 58,032.19 was common to all variant models. Using these initial data and shape models, we propagated all 687 clone realizations of 2012 TC4 backward in time to the epoch MJD 56,209.88, characteristic of the 2012 October observations. We used our numerical approach described above with both gravitational and radiation (YORP) torques included. For the latter, we assumed an effective size D = 10 m (corresponding to a sphere of the same volume of the models) and bulk density ρ = 1.4 g cm⁻³. The abovementioned rescaling rule, namely, invariance to ρ D², may allow us to transform our results to other combinations of D and ρ values.

The evolutionary tracks of angular momentum L and energy E secular changes due to the YORP torques mostly resemble those from Figure 7. For some model variants, which were more different from the nominal one, the slope of the overall secular change in L and/or E was shallower or steeper. At the initial and final epochs of our simulations, we determined P_ϕ and P_ψ periods from a short numerical simulation of a free-top model. We also verified that the P_ψ values correspond exactly to those provided by the analytical formula (Equation (A5)). Figure 8 shows our results. The blue histograms are for the initial data, i.e., the 2017 October rotation state. Because the observations were numerous and of good quality, the model variants differ only slightly, and the P_ϕ and P_ψ distributions are tight (they also match those from Figure 4). The red histograms in Figure 8 were determined from the last epoch of our numerical runs and correspond to the predicted rotation state of 2012 TC4 in 2012 October. We note that both periods appreciably changed, as already anticipated from the preliminary simulation shown in Figure 7. These P_ϕ and P_ψ distributions are obviously less tight than the initial ones, reflecting the different evolutionary tracks of the individual clone variants. Our prime interest is to compare the red distribution in Figure 8 (from simulations) to the red distributions in Figure 4 (from the 2012 observations; also highlighted with a dashed line for the P_ϕ period). First, we note that the match of the P_ϕ periods is surprisingly good. The mean value of 8.495 minutes of the observations corresponds rather well to the mean value of 8.497 minutes of the simulated data, which have a comfortably large dispersion of 0.003 minutes to overlap with the observed data; in fact, the shift in P_ϕ is even larger than required. Interestingly, the comparison is not as good in the P_ψ period. The model-predicted value of 27.59 ± 0.02 minutes (formal uncertainty) is short to explain the observations, which provide, on average, 27.87 minutes. Still, the model indicates a significant shift from the initial value of 27.5070 ± 0.0002 minutes. Nevertheless, to reach the value from the 2012 observations, the shift would need to be about 3.7 times larger. We do not know the reason for this difference; in all likelihood, it is related to the misbehavior in the rotational energy evolution (see above). We suspect that the overly simple modeling of the YORP effect, such as the surface thermal inertia and/or the unresolved small-scale shape irregularities, plays an important role here (note that a factor of 3 between the observations and model prediction was also seen in the cases where YORP was detected for asteroids rotating about the principal axis of the inertia tensor). The fact that some deeper aspects of the model are not characterized well, as witnessed by the opposite sign of the energy evolution, implies that our results cannot be easily reconciled with the observations by a simple rescaling of size D and bulk density ρ. We have verified that the accumulated shifts in both the P_ϕ and P_ψ periods are proportional to ρ D². Thus, the P_ψ mismatch could be explained by assuming that 2012 TC4's size is ≃5.2 m, but this would produce a factor of ≃3.7 inconsistency in the P_ϕ period (making the modeled value larger than observed). Instead, we believe that some missing details of the YORP modeling, which result in different effects on P_ϕ and P_ψ, are responsible for the difference.

Zoom In Zoom Out Reset image size

4.3.2. Model Based on a Combination of 2012 and 2017 Data

For the sake of comparison, we also repeated our analysis using models based on a combination of the observations in 2012 and 2017 (Section 3.3). Obviously, at each of these epochs, we considered different parameters of the rotation state, but now we enforce that the same shape model is used for both data sets. This solution gives us an opportunity to consider two sets of initial conditions for our simulation, in both 2012 and 2017, and analyzes the predictions in the complementary epoch (integrating the rotation model once backward in time and once forward in time). Obviously, in all cases, our model includes the gravitational and radiation torques, D = 10 m effective size, and ρ = 1.4 g cm⁻³ bulk density, as before (needed for the evaluation of the radiation torques).

We start with the case of the initial data in 2017 October and backward-in-time integration. This is directly comparable with the results above, when only observations in 2017 were used. However, the models are slightly different in all aspects (initial conditions and shape) because now the 2012 observations play a role in their construction. The results are shown in Figure 9. While slightly different than in Figure 8, the principal outcome is the same: (i) a fairly satisfactory prediction for the P_ϕ period and (ii) a too-small change in the P_ψ period.

Zoom In Zoom Out Reset image size

We next consider the opposite case, namely, the initial condition in 2012 October and model propagation forward in time to 2017 October. The results are shown in Figure 10. Obviously, here the red histograms (from the initial data in 2012 October) are more constrained than the blue histograms (resulting from rotation state vectors propagated using our model to 2017 October). The latter are more dispersed than the red histograms in Figure 9 because the less numerous and accurate observations in 2012 constrain the models with lower accuracy. Nevertheless, the principal features of the solution are still present: (i) the P_ϕ period changed adequately for the majority of cases, while (ii) the P_ψ period changed too little.

Zoom In Zoom Out Reset image size

4.3.3. In Summary

The individual solutions of 2012 TC4's rotation parameters in 2012 and 2017 indicate that periods P_ψ and P_ϕ decreased in the latter epoch (Figure 4). In terms of an osculating approximation with a free-top model, it implies that the wobbling motion of the angular momentum vector L in the body-fixed frame moved toward the fundamental mode of its direction along the shortest axis of the inertia tensor. Note that the trajectory of L in the body-fixed frame is uniquely parameterized with p = 2BE/L² (see, e.g., Appendix A and Landau & Lifshitz 1969). In quantitative terms, the change from the 2012 to 2017 state is expressed by δp ≃ − 6.2 × 10⁻³ (Table 2). The average rate over a δt ≃ 5 yr interval would thus be δp/δt ≃ −3.9 × 10⁻¹¹ s⁻¹.

In the YORP model presented above, the change in p was a composition of changes in both the energy E and angular momentum L. In fact, both E and L increased from 2012 to 2017 (Figure 7), but the composite effect was a decrease in p, in this case, δp ≃ − 4.6 × 10⁻³, a similar value to that directly determined from osculating L and E above. Other processes may lead to approximately the same results by producing a different combination of energy and angular momentum changes. For instance, the effects of material inelasticity directly result in energy dissipation while preserving angular momentum. In this case, both E and p decrease with the direct relation δE ≃ (L²/2B)δp.

In order to explore whether the observed effect of 2012 TC4's spin change could even be plausibly matched by internal energy dissipation, we used the model presented by Breiter et al. (2012). These authors assumed a fully triaxial geometry of the body but restricted their analysis of energy dissipation to the empirical description with a quality factor Q (see an alternative model of Frouard & Efroimsky 2017, where the authors described the energy dissipation using a Maxwell viscous liquid but allowed only a biaxial shape of the body). With these assumptions, they expressed the secular (i.e., wobbling cycle–averaged) rate of energy change in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{\delta E}{\delta t}\simeq -\displaystyle \frac{{a}^{4}\rho m{{\rm{\Omega }}}^{5}}{\mu Q}\,{\rm{\Psi }},\end{eqnarray} \tag{ 1 }$$

where a is the semimajor axis of the body's triaxial approximation, ρ is its density, m is its mass, Ω = L/C, and Ψ is a rather complicated factor depending on the nutation angle θ_nut (i.e., the tilt between L and the shortest body axis; we find that θ_nut oscillates between ≃16° and ≃46°, with a mean of ≃30°), body axis ratios, and Lamé coefficients (see Breiter et al. 2012, Section 4.3). Finally, μ is the Lamé shear modulus (rigidity) and Q is the quality factor, empirically expressing the energy dissipation per wobbling cycle. The product μQ has been characteristic of studies involving energy dissipation in planetary science since the pioneering work of Burns & Safronov (1973; see, however, Prendergast 1958). While highly uncertain, the typical values of this parameter for asteroids range between 10¹¹ and 5 × 10¹² Pa (e.g., Harris 1994). Using δE/δt ≃ (L²/2B)δp/δt with the abovementioned δp/δt value, we can now use Equation (1) to infer what values of μQ would be needed to explain the change in 2012 TC4's tumbling state in between the 2012 and 2017 close approaches. The remaining unknown parameter is Ψ, which we estimate to be ≃(1−5) × 10⁻³. Plugging this value into Equation (1), we find μQ ranging from 4 × 10⁵ to 4 × 10⁶ Pa. Remarkably, such values are 4–5 orders of magnitude smaller than the usually adopted estimates. Therefore, unless the energy dissipation by internal friction is extraordinarily high (and thus the μQ value is very small), this process cannot explain the observed rotation change of 2012 TC4. Note, additionally, that we considered the derivation of the needed μQ value within the energy dissipation model as a useful exercise to match the energy change. Such a model, however, would not explain the observed angular momentum change.

Another alternative process to the radiation torques is that of an impact by an interplanetary meteoroid. However, in the following, we provide an argument that the likelihood of this happening at the level needed to explain the 2012 TC4 data is again very small. To that end, we used information from Bottke et al. (2020). They determined the meteoroid flux on a small asteroid, (101955) Bennu, using the state-of-the-art model MEM-3, allowing them to predict the parameters of meteoroids impacting a target body orbiting between Mercury and the asteroid belt (e.g., Moorhead et al. 2020). Note that the orbit of Bennu is similar to 2012 TC4, and we shall neglect the small flux differences that could result from a small orbital dissimilarity of these two objects (if anything, the flux would be slightly larger on Bennu because of its closer orbit to the Sun). In their Figure 2, Bottke et al. (2020) showed that 5 mg interplanetary particles, approximately 2 mm in size, strike Bennu with a frequency of ≃60 yr^–1 with a median impact velocity of a little less than ≃30 km s⁻¹. For 2012 TC4, we only need to rescale this number to account for a much smaller size; Bennu is an ≃500 m asteroid, while the characteristic size of 2012 TC4 is only ≃10 m. Therefore, the 5 mg flux on 2012 TC4 is about ≃2.4 × 10⁻² yr^–1. The chances of being hit by such a particle in 5 yr is therefore merely ≃0.12. However, even if it had happened, the dynamical effect would be minimum. Estimating the change in rotational angular momentum L plainly by δL/L ≃ (m/M)(v_imp/Rω), where m and M are the masses of the particle and 2012 TC4, v_imp is the impact velocity, and R and ω are 2012 TC4's characteristic radius and rotational frequency, we would obtain δL/L ≃ 2 × 10⁻⁷. At least a 10,000 times larger effect would be needed to approach the level observed for 2012 TC4, and this would require an impacting particle at least 20 times larger, i.e., 5 cm or more. Because MEM-3 incorporates the flux dependence on mass from Grün et al. (1985), thus ∝ m^−4/3, the chances that 2012 TC4 was hit by a 5 cm meteoroid between 2012 and 2017 is only ≃6× 10⁻⁷. We thus conclude that the chances that the observed effect was produced by an impact of an interplanetary particle are negligibly small and may be discarded.