The Changing Rotational Light-curve Amplitude of Varuna and Evidence for a Close-in Satellite

Trans-Neptunian objects (TNOs) are solar system bodies that orbit the Sun with larger semimajor axes than that of Neptune, formed quite outside the so-called "snow line" where the temperature of the protoplanetary disk was low enough to allow the survival of molecules of chemical compounds with low sublimation points. Due to the ample distances that separate TNOs from the Sun, they have suffered less chemical processes than other solar system bodies and may preserve material from the nebula that generated the solar system. Hence, TNOs yield important information about the solar system formation and its evolution.

Currently, observational searches have found ∼2500 TNOs (according to the data from MPC⁷ ), although dynamical models suggest that the TNO population located between 30 and 50 au from the Sun is around 10⁵ objects with diameters larger than 100 km (Trujillo et al. 2001; Petit et al. 2008, 2011).

20000 Varuna is one of the most interesting TNOs due to its peculiar physical properties. It rotates relatively fast with a period of 6.34358 ± 0.00005 hr (Belskaya et al. 2006; Santos-Sanz et al. 2013), producing a double-peaked rotational light curve dominated by its shape (e.g., Jewitt & Sheppard 2002; Lellouch et al. 2002; Hicks et al. 2005; Belskaya et al. 2006). The rotational light-curve amplitude reported in these studies is quite large, ∼0.45 mag, indicating that the body is highly elongated (i.e., b/a < 0.66, being a and b the larger semiaxes of a triaxial body). The fast period and the elongated shape require a density of ∼1000 kg m⁻³ (assuming hydrostatic equilibrium; Chandrasekhar 1987). This is somewhat high compared to other TNOs of similar sizes, considering that Varuna's equivalent diameter is ∼700 km, given by thermo-physical models using Herschel Space Telescope measurements (Lellouch et al. 2013); see the supplementary material in Ortiz et al. (2012) to compare the density of objects with similar sizes. However, a stellar occultation by Varuna, detected in 2010, results in a long chord of 1004 km (Sicardy et al. 2010). This value is somewhat in tension with the one given by Lellouch et al. (2013), but fits better with Varuna's density.

Motivated by all of this, we have been monitoring Varuna for nearly two decades now, resulting in a collection of data with several rotational light curves (which we will simply call light curves in the following). In this work, we present new observations of Varuna at different epochs (Section 2). The results from these observations, together with other observations from the literature, are presented in Section 3. Section 4 contains the analysis of the set of light curves in order to obtain Varuna's pole orientation and shape. Section 5 displays the evidence for a possible satellite. A brief summary is presented in Section 6.

We carried out eight observing campaigns from 2005 to 2019 using four different telescopes. A log of observations is shown in Table 1.

The 1.5 m telescope at the Sierra Nevada Observatory in Granada (Spain) has a CCD VersArray with 2048 × 2048 pixels, while the field of view (FoV) is 792 × 792 and the image scale is 0232 pixel⁻¹. Observations were obtained in 2 × 2 binning mode. The typical seeing was ∼15, being the point-spread function well sampled for the photometric goals.

The 1.23 m and the 2.2 m telescopes at Calar Alto Observatory are located in Almería (Spain). The 1.23 m telescope uses a DLR-MKIII camera with 4096 × 4096 pixels, while the FoV and the image scale are 215 × 215 and 0315 pixel⁻¹, respectively. At the 2.2 m telescope we used two different instruments: the Bonn University Simultaneous CAmera (BUSCA) and the Calar Alto Faint Object Spectrograph (CAFOS). BUSCA allows the simultaneous direct imaging of the same sky area in four colors; each CCD has 4096 × 4096 pixels, with an FoV of 12' × 12' and an image scale of 0176 pixel⁻¹. We used 2 × 2 binning mode. CAFOS is equipped with a 2048 × 2048 pixel CCD; we used the SITe#1d chip, which produces a 053 pixel⁻¹ image scale and a circular FoV of 16'.

The Telescopio Nazionale Galileo (TNG), with 3.58 m diameter, is located at the Roque de los Muchachos Observatory in Canary Islands (Spain). We used the Device Optimized for the LOw RESolution (DOLORES) with a detector of 2048 × 2048 pixels. The image scale is 0252 pixel⁻¹, which yields an FoV of 86 × 86.

Bias and twilight flat-field calibration images were taken at the beginning of each observation night. The science images were processed by subtracting a median bias and dividing by a median flat-field. Using the Interactive Data Language, we developed our own code to perform these reduction steps, as well as the synthetic aperture photometry for the time series data.

We performed the synthetic aperture photometry using different aperture radii, as described in Fernández-Valenzuela et al. (2016), and choosing the one in which the dispersion of the photometry is minimal. To construct the light curves, Varuna was compared to a set of reference stars located as close as possible to the object. Whenever possible during each campaign, we employed the same set of comparison stars each night in order to minimize systematic photometric errors.

From the time series observations obtained between 2005 and 2019, we have built eight different light curves. Additionally, we used data published in the literature by Jewitt & Sheppard (2002), Lellouch et al. (2002), and Hicks et al. (2005), thus incorporating three additional light curves to our study from previous years (Figure 1). All the light curves were corrected from light travel time. Because Varuna's body is assumed to have an ellipsoidal shape (e.g., Jewitt & Sheppard 2002; Lellouch et al. 2002), data from each light curve were fitted to a Fourier series $m={{\rm{\Sigma }}}_{i}\left[{a}_{i}\sin (2i\pi \phi )+{b}_{i}\cos (2i\pi \phi )\right]$, where m is the theoretical value of the relative magnitude obtained from the fit, ϕ is the rotational phase (calculated as the fractional part of (JD − JD₀)/P, where JD is the Julian Date, JD₀ = 2451957.0 is the initial Julian Date, and P is the rotation period in days), and (a_i, b_i) are the coefficients of the Fourier function (with i = 0, 1, 2, ...). In our specific case, we used up to second-order (i = 2) or up to fourth-order (i = 4) Fourier functions. The second order is the minimum order that allows a double-peaked fit; however, higher orders take into account small deviations on inhomogeneous objects and can be used to fit light curves that are highly sampled. The 2001, 2002A, 2011, 2018, and 2019 light curves were fitted to a fourth-order function, while the remaining light curves were fitted to a second-order Fourier function, because the number of data points in those runs was not large enough to use a higher order. Data were folded using Varuna's rotation period of 6.343572 ± 0.000006 hr, which is obtained using the Lomb periodogram analysis of all our data, in agreement with the previous one reported in Belskaya et al. (2006). The peak-to-valley amplitudes Δm (amplitudes in the following) of each light curve are given by the absolute maximum and minimum produced by the fits.

Zoom In Zoom Out Reset image size

Table 1 contains the results from the fit to each light curve, i.e., the amplitude and the dispersion of the residuals of the Fourier function fit to the observational data. One of the evident results is that the amplitude has changed considerably along these 19 years, with an increase of ∼0.13 mag. Table 2 presents all the relative photometry observations from 2005 to 2019.

The aspect angle of a solar system body, formed by the line of sight and the rotational axis of the body, changes due to the orbital motion of the body around the Sun. This change produces a variation in the geometric cross section of the object's rotational modulation that is translated into a variation in the amplitude which, for a triaxial ellipsoid, is described as

$$\begin{eqnarray}&&{\rm{\Delta }}m=-2.5\mathrm{log}\left[\displaystyle \frac{b}{a}{\left(\displaystyle \frac{{a}^{2}{\cos }^{2}(\delta )+{c}^{2}{\sin }^{2}(\delta )}{{b}^{2}{\cos }^{2}(\delta )+{c}^{2}{\sin }^{2}(\delta )}\right)}^{1/2}\right],\end{eqnarray} \tag{ 1 }$$

where δ is the aspect angle and a, b, c are the semiaxes of a triaxial body (with a > b > c; Tegler et al. 2005). The aspect angle of the body can be written as

$$\begin{eqnarray}\begin{array}{rcl}\delta & = & \displaystyle \frac{\pi }{2}-\arcsin \left[\sin ({\beta }_{{\rm{e}}})\sin ({\beta }_{{\rm{p}}})\right.\\ & & +\left.\,\cos ({\beta }_{{\rm{e}}})\cos ({\beta }_{{\rm{p}}})\cos ({\lambda }_{{\rm{e}}}-{\lambda }_{{\rm{p}}})\right],\end{array}\end{eqnarray} \tag{ 2 }$$

where λ_e and β_e are the ecliptic longitude and latitude of Varuna's-centered reference frame (given by the ephemeris), and λ_p and β_p are the ecliptic longitude and latitude of the pole orientation (Schroll et al. 1976). These coordinates change with the observing epochs. In order to fit the observational data to Equation (1), we carried out a grid search for the free parameters in the equations: b/a, c/b, λ_p, and β_p, which gave theoretical values for Δm (as done in Fernández-Valenzuela et al. 2017). In a first step, we assumed hydrostatic equilibrium, thus reducing the number of free parameters to three, as c/b is then related to b/a (Chandrasekhar 1987). We explored b/a in the range [0.44, 0.60] at intervals of 0.02. The upper limit of 0.60 is imposed by the measurement of the largest observed amplitude up to date (see Table 1). The lower limit of 0.44 is imposed by theoretical models of rotational equilibrium configurations of strengthless bodies (Leone et al. 1984). We explored the possible values of λ_p and β_p on the entire sky at intervals of 0.5°.

The goodness of the fit is evaluated using the ${\chi }_{\mathrm{pdf}}^{2}$ test. The possible solutions are those that have values of ${\chi }_{\mathrm{pdf}}^{2}$ within the range ${\chi }_{\mathrm{pdf},\min }^{2}+0.584$ (from the χ² distribution with a 0.9 level of confidence for 3 degrees of freedom; Feller 1971). As a result, the subsequent range of values for the axis ratio b/a is [0.56, 0.60], where the lower limit is given by the goodness of the fit. Considering the results from the different axis ratios, we constrained the pole orientation with ${\lambda }_{{\rm{p}}}\in [43.0,83.5]^\circ$ and ${\beta }_{{\rm{p}}}\in [-70.0,-58.0]^\circ$ and ${\lambda }_{{\rm{p}}}\in [209,220]^\circ$ and ${\beta }_{{\rm{p}}}\in [-49.0,-35.5]^\circ$ (with their complementary directions also possible).

Specifically, we obtained the ${\chi }_{\mathrm{pdf},\min }^{2}=1.904$ for the following parameters: semiaxis ratio b/a = 0.60 (with c/b = 0.72 imposed by the hydrostatic equilibrium), and pole orientation λ_p = 53° ± 10° and β_p = −64° ± 6° (the complementary direction λ_p = 233° and β_p = 64° is also possible for the same ${\chi }_{\mathrm{pdf}}^{2}$ value). The errors of λ_p and β_p have been obtained considering the values that fit within the interval of ${\chi }_{\mathrm{pdf},\min }^{2}+0.584$ for this specific axis ratio. Figure 2 shows the best fit given by the above values (left panel) and a ${\chi }_{\mathrm{pdf}}^{2}$ map as a function of λ_p and β_p (right panel) for the specific value of b/a = 0.60. The solution in the range of ${\lambda }_{{\rm{p}}}\in [209,211]^\circ$ and ${\beta }_{{\rm{p}}}\in [-42.0,-35.5]^\circ$ can be discarded making use of results from the stellar occultation that occurred in 2010 (Sicardy et al. 2010). Even if a detailed analysis will be published elsewhere, we can advance that this yields restricted intervals for the values of λ_p and β_p. It turns out that the first possibility ${\lambda }_{{\rm{p}}}=53^\circ \pm 10^\circ$, β_p = −64° ± 6° (as well as its complementary direction) lies inside these intervals, for which the resulting area of Varuna does not produce a geometric albedo being unnaturally small,⁸ but this is not the case for ${\lambda }_{{\rm{p}}}\in [209,211]^\circ$, ${\beta }_{{\rm{p}}}\in [-42.0,-35.5]^\circ$, for which the size of Varuna would be too large, corresponding to an albedo well below the lower limit of 0.04 for TNOs.

Zoom In Zoom Out Reset image size

In a second step, we relaxed the assumption of hydrostatic equilibrium, adding the parameter c/b to the grid search in the range [0.5, 1.0] using steps with a length of 0.1. However, this does not lead to appreciable differences. The smallest value of ${\chi }_{\mathrm{pdf}}^{2}$ found is 1.905, which results in values of λ_p = 54° and β_p = −65°, with the axis ratios b/a = 0.60 and c/b = 0.70, which is the equivalent result to the case of hydrostatic equilibrium. Currently, from the best fit, Varuna is close to an "edge-on" configuration (see the left panel in Figure 2), where the aspect angle is close to 90°. Therefore, the c-axis, which is the one that would determine whether Varuna is in hydrostatic equilibrium or not, does not produce a strong effect on the light-curve shape, being difficult to distinguish between both scenarios.

There also exists the possibility that Varuna could be a contact binary, as happened with 2014 MU₁₉ (Stern et al. 2019). However, this scenario was explored by Jewitt & Sheppard (2002) and was considered very unlikely for Varuna. Note that contact binaries have specific V-shaped light curves (Thirouin & Sheppard 2017; Thirouin et al. 2017), which is not the case of Varuna.

Even though the signal-to-noise ratio (S/N) of the data in most of the runs was good (well above 50), the residuals of the Fourier fits were considerably scattered, indicating that the model may have been missing some aspect. We focused on the data obtained during 2019, with the best S/N (∼80) for Varuna. In Figure 3, we plot the relative magnitude of each day versus Julian date. While the four light curves present a good fit, at localized points there are deviations of the model with respect to the data. Moreover, when merging the four days, the quality of the assembled light curve decreases considerably, as the scatter of residuals increases with respect to single night fits and the parameters of the fits change slightly from night to night. A satellite might be responsible for changes in the shape of the main Varuna's light curve. In order to test the potential presence of a rotating satellite that could be generating an additional periodicity in the photometry, we used the residuals from the fit of the 2011, 2014A, 2014B, 2018, and 2019 light curves (which are the ones with higher quality that are better sampled) and searched for a period using the Lomb periodogram technique (Lomb 1976) as implemented in Press et al. (1992). This technique has been largely applied to detect satellites in asteroids (e.g., Pravec & Hahn 1997; Pravec et al. 2002; Margot et al. 2015). In Figure 4 we plotted the Lomb periodogram of the residuals (left panel), which gives the maximum spectral power of 45 for the frequency of 2.003 cycles/day (11.9819 hr), and the residuals folded to that period (right panel), where the blue line shows a one-order Fourier function fit. The amplitude obtained from this fit, 0.04 mag, is related to the rotational modulation of the satellite. Assuming that the satellite's light curve is due to an albedo spot and that its spin and orbit are synchronized (we will check below that this is a reasonable assumption), the distance d at which the satellite orbits from Varuna is given by the equation:

$$\begin{eqnarray}&&{d}^{3}=\displaystyle \frac{{T}_{{\rm{S}}}^{2}{GM}}{4\pi },\end{eqnarray} \tag{ 3 }$$

where T_S is the satellite's orbital period (in seconds), G is the gravitational constant and M = M_V + M_S is the mass of the system. Varuna's mass, ∼10²¹ kg, is estimated taking into account a density of 1100 kg m⁻³, assuming hydrostatic equilibrium in a Jacobi ellipsoid with axis ratio $b/a\in [0.60-0.56]$, and volume-equivalent diameter of 700 km. This results in an estimation of d ∼ 1300 km, which places the satellite outside the Roche limit,

$$\begin{eqnarray}&&{d}_{{\rm{R}}}={a}_{{\rm{V}}}{\left(2\displaystyle \frac{{\rho }_{{\rm{V}}}}{{\rho }_{{\rm{S}}}}\right)}^{1/3},\end{eqnarray} \tag{ 4 }$$

where a_V ≃ 550 km and ρ_V are Varuna's largest semiaxis and density, respectively, and ρ_S is the satellite's density. Assuming that the satellite has a density of ∼300 kg m⁻³ (similar to 2014 MU₆₉'s density; Stern et al. 2019), this results in a Roche limit of 1000 km.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We should note that the satellite could have an elongated shape with a double-peaked light curve. This implies that the rotation period of the satellite would be double (∼24 hr), and so would its orbital period. This would locate the satellite ∼2000 km from Varuna, with no strong implications for the following calculations. The angular distance between Varuna and the satellite is then smaller than the resolution of the Wide Field Camera 3 installed at the Hubble Space Telescope, thus making it impossible to resolve the system.

The time required by the satellite to be tidally locked, t_lock, can be obtained using the equation given by Hubbard (1984),

$$\begin{eqnarray}&&{t}_{\mathrm{lock}}=\displaystyle \frac{2\pi {M}_{{\rm{V}}}{d}^{6}}{3{k}_{{\rm{V}}}{{GM}}_{{\rm{S}}}^{2}{a}_{{\rm{V}}}^{3}{T}_{{\rm{V}}0}\delta },\end{eqnarray} \tag{ 5 }$$

where M_V is Varuna's mass, k_V is the secular Lomb number, which has a value of 3/2 for homogeneous bodies, M_S is the mass of the satellite, T_V0 is the initial rotation rate of Varuna, and δ is expressed as arctan(1/Q) with Q is the dissipation. Assuming a dissipation of Q = 100, and that T_V0 was greater than or equal to T_V currently measured, it is straightforward to calculate the dimensionless number t_lock/t_ss with t_ss ≃ 1.5 × 10¹⁷ s the age of the solar system, obtaining ${t}_{\mathrm{lock}}/{t}_{\mathrm{ss}}\lesssim 5\times {10}^{9}\,{\rm{kg}}\times \left({M}_{{\rm{V}}}/{M}_{{\rm{S}}}^{2}\right)$. Using Varuna's largest semiaxis and its density, it follows that, as long as M_S/M_V ≫ 10⁻⁶, then the dimensionless number t_lock/t_ss, is small. This validates the assumption that the satellite is tidally locked, as the time required by the the satellite to enter this configuration is much smaller than the age of the solar system.

The value obtained for the pole orientation (Section 4) places the body almost in an "edge-on" configuration, implying that mutual events produced by the satellite may occur in the coming years, assuming that the satellite's orbit is in Varuna's equatorial plane.

Time-series photometry of Varuna at different epochs have revealed a remarkable change in the amplitude of the light curve of ∼0.13 mag. A geometric model that explains this change as the rotational variability due to Varuna's shape provides a range of solutions for its pole orientation of λ_p ∈ [43.0, 83.5]° and ${\beta }_{{\rm{p}}}\in [-70.0,-58.0]^\circ$ and its axes ratio b/a ∈ [0.56, 0.60]. Nowadays, only around four TNOs have well determined pole orientations (Pluto, Charon, Haumea, and 2014 MU₆₉), with three of them visited by spacecraft. Additionally, we noticed that Varuna's light curves have deviations from the expected behavior of a simple ellipsoidal body, and we analyzed the periodogram of the residuals of the photometric fits to the data, finding a peak of high spectral power for a period of 11.9819 hr. This appears consistent with the idea of the existence of a close-in satellite, separated from Varuna a distance of ∼1300 km (depending on Varuna's mass, which is not accurately known). Until now, satellite discoveries in the trans-Neptunian region have been only accomplished using direct imaging (Noll et al. 2008). With this work, we have illustrated that photometric techniques can be used to reveal the existence of a satellite for an object in the trans-Neptunian region. The light-curve technique has the potential to reveal objects much closer-in than the direct imaging technique, helping to remove biases in our current knowledge of the population of binaries or bodies with satellites.

E.F.V. acknowledges UCF 2017 Preeminent Postdoctoral Program (P³). Part of the research leading to these results has received funding from the European Unions Horizon 2020 Research and Innovation Programme, under Grant Agreement No. 687378 (SBNAF). We would like to acknowledge financial support by the Spanish grant AYA-2017-84637-R and the financial support from the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV- 2017-0709). This work is partially based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. We thank the director of the Telescopio Nazionale Galileo for allocation of Director's Discretionary Time. We acknowledge Raúl Carballo-Rubio for providing useful comments that helped improve this manuscript and Pedro Gutiérrez for helpful discussions.