The changing rotational light-curve amplitude of Varuna and evidence for a close-in satellite

From CCD observations carried out with different telescopes, we present short-term photometric measurements of the large trans-Neptunian object Varuna in 10 epochs, spanning around 19 years. We observe that the amplitude of the rotational light-curve has changed considerably during this period of time from 0.41 to 0.55 mag. In order to explain this variation, we constructed a model in which Varuna has a simple triaxial shape, assuming that the main effect comes from the change of the aspect angle as seen from Earth, due to Varuna's orbital motion in the 19-year time span. The best fits to the data correspond to a family of solutions with axial ratios b/a between 0.56 and 0.60. This constrains the pole orientation in two different ranges of solutions presented here as maps. Apart from the remarkable variation of the amplitude, we have detected changes in the overall shape of the rotational light-curve over shorter time scales. After the analysis of the periodogram of the residuals to a 6.343572 h double-peaked rotational light-curve fit, we find a clear additional periodicity. We propose that these changes in the rotational light-curve shape are due to a large and close-in satellite whose rotation induces the additional periodicity. The peak-to-valley amplitude of this oscillation is in the order of 0.04 mag. We estimate that the satellite orbits Varuna with a period of 11.9819 h (or 23.9638 h), assuming that the satellite is tidally locked, at a distance of ~ 1300 km (or ~ 2000 km) from Varuna, outside the Roche limit.


INTRODUCTION
Trans-Neptunian objects (TNOs) are solar system bodies that orbit the Sun with larger semi-major axis than that of Neptune, formed quite outside the so-called "snow line" where the temperature of the proto-planetary 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 1 ), although dynamical models suggest that the TNO population located between 30 − 50 au from the Sun is around 10 5 objects with diameters larger than 100 km (Trujillo et al. 2001;Petit et al. 2008Petit et al. , 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 h (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 semi-axes of a triaxial body). The fast period and the elongated shape require a density of ∼ 1000 kg m −3 (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 thermophysical models using Herschel Space Telescope measurements (Lellouch et al. 2013); see the supplementary material in Ortiz et al. (2012) to compare the density of ob-jects 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 (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.

OBSERVATIONS AND DATA REDUCTION
We carried out eight observing campaigns from 2005 to 2019 using four different telescopes. A log of observations is shown in online table 1.
The 1.5-m telescope at Sierra Nevada Observatory in Granada (Spain) has a CCD VersArray with 2048 × 2048 pixels, while the field of view (FoV) is 7.92 ′ × 7.92 ′ and the image scale is 0.232 ′′ /pixel. Observations were obtained in 2 × 2 binning mode. The typical seeing was ∼ 1.5 ′′ , 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 21.5 ′ × 21.5 ′ and 0.315 ′′ /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 a FoV of 12 ′ × 12 ′ and an image scale of 0.176 ′′ /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 0.53 ′′ /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 0.252 ′′ /pixel which yields a FoV of 8.6 ′ × 8.6 ′ .
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 (IDL), 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.

RESULTS FROM OBSERVATIONS
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 = Σ i [a i sin(2iπφ) + b i cos(2iπφ)], 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 0 )/P , where JD is the Julian Date, JD 0 = 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 secondorder (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 2001The , 2002A, 2011The , 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 h, 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. 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. Online table 2 presents all the relative photometry observations from 2005 to 2019.

VARUNA'S POLE ORIENTATION AND SHAPE
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  (2002); Lellouch et al. (2002) and Hicks et al. (2005), respectively. The lines over plotted to the data points represent the fit of the observational data to the Fourier series. Error bars are not shown to avoid a cluttered plot; the data with the uncertainties are given in the online table 2.
a triaxial ellipsoid, is described as ∆m = −2.5 log b a a 2 cos 2 (δ) + c 2 sin 2 (δ) b 2 cos 2 (δ) + c 2 sin 2 (δ) where δ is the aspect angle and a, b, c are the semi-axes of a triaxial body (with a > b > c; Tegler et al. 2005). The aspect angle of the body can be written as 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 the 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  Lellouch et al. (2002) c Data from Hicks et al. (2005) 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 χ 2 pdf test. The possible solutions are those that have values of χ 2 pdf within the range χ 2 pdf,min + 0.584 (from the χ 2 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 λ p ∈ [43.0, 83.5] • and β p ∈ [−70.0, −58.0] • and λ p ∈ [209, 220] • and β p ∈ [−49.0, −35.5] • (with their complementary directions also possible).
Specifically, we obtained the χ 2 pdf,min = 1.904 for the following parameters: semi-axis 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 χ 2 pdf value). The errors of λ p and β p have been obtained considering the values that fit within the interval of χ 2 pdf,min + 0.584 for this specific axis ratio. Figure 2 shows the best fit given by the above values (left panel) and a χ 2 pdf 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 λ p ∈ [209, 211] • and β p ∈ [−42.0, −35.5] • 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 λ p = 53 • ± 10 • , β 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 2 , but this is not the case for λ p ∈ [209, 211] • , β p ∈ [−42.0, −35.5] • , for which the size of Varuna would be too large, corresponding to an albedo well below the lower limit of 0.04 for TNOs.
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 χ 2 pdf 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 "edgeon" configuration (see 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 19 (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 et al. 2017;Thirouin and Sheppard 2017), which is not the case of Varuna.

INDICATIONS FOR VARUNA'S CLOSE-IN SATELLITE
Even though the signal-to-noise ratio (SNR) 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 SNR (∼ 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 responsable 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 and 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 and 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 h), 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: 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 21 kg, is estimated taking into account a density of 1100 kg m −3 , assuming hydrostatic equilibrium in a Jacobi ellipsoid with axis ratio b/a ∈ [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, where a V ≃ 550 km and ρ V are Varuna's largest semi-axis and density, respectively, and ρ S is the satellite's density. Assuming that the satellite has a density of ∼ 300 kg m −3 (similar to 2014 MU 69 's density; Stern et al. 2019), this results in a Roche limit of 1000 km. We should note that the satellite could have an elongated shape with a doublepeaked light-curve. This implies that the rotation period of the satellite would be double (∼ 24 h), and so 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 imposible 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), 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 or equal than T V currently measured, it is straightforward to calculate the dimensionless number t lock /t ss with t ss ≃ 1.5 × 10 17 s the age of the solar system, obtaining t lock /t ss 5 × 10 9 kg × (M V /M 2 S ). Using Varuna's largest semi-axis and its density, it follows that, as long as M S /M V ≫ 10 −6 , 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.  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.

ACKNOWLEGEMENTS
EFV acknowledges UCF 2017 Preeminent Postdoctoral Program (P 3 ). 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 improving this manuscript and Pedro Gutiérrez for helpful discussions.