A multi-chord stellar occultation by the large trans-Neptunian object (174567) Varda

We present results from the first recorded stellar occultation by the large trans-Neptunian object (174567) Varda that was observed on September 10$^{\rm th}$, 2018. Varda belongs to the high-inclination dynamically excited population, and has a satellite, Ilmar\"e, which is half the size of Varda. We determine the size and albedo of Varda and constrain its 3D shape and density. Thirteen different sites in the USA monitored the event, five of which detected an occultation by the main body. A best-fitting ellipse to the occultation chords provides the instantaneous limb of the body, from which the geometric albedo is computed. The size and shape of Varda are evaluated, and its bulk density is constrained, using Varda's mass known from previous works. The best-fitting elliptical limb has semi-major (equatorial) axis of $(381 \pm 3)$km and an apparent oblateness $0.043\pm0.036$ corresponding to an apparent area-equivalent radius $R'_{\rm equiv}= (373\pm8)$km and geometric albedo $p_v=0.097\pm 0.004$ assuming a visual absolute magnitude $H_V=3.81\pm0.01$. Using three possible rotational periods for the body (4.76~h, 5.91~h, and 7.87~h), we derive corresponding MacLaurin solutions. Furthermore, given the low-amplitude ($0.06\pm0.01$) mag of the single-peaked rotational light-curve for the aforementioned periods, we consider the double periods. For the 5.91~h period (the most probable) and its double (11.82~h), we find bulk densities and true oblateness of $\rho=(1.52\pm0.05)$ g cm$^{-3}$, $\epsilon=0.232\pm0.036$ and $\rho=(1.25\pm0.04)$ g cm$^{-3}$, $\epsilon=0.079\pm0.044$. However, it must be noted that the other solutions cannot be excluded just yet.


Introduction
Stellar occultations are one of the most accurate ground-based methods to directly determine the size and shape of solar system objects, down to kilometric accuracies. In addition, these rare events can be used, amongst others, to reveal the existence of atmospheres, cometary activity, satellites or ring systems around minor planets (Ortiz et al. 2020). We can cite here the discovery of rings around the Centaur (10199) Chariklo (Braga-Ribas et al. 2014) and the dwarf-planet (136108) Haumea ) as well as the detection of (87) Sylvia's satellites (Vachier et al. 2019).
In this paper, we report results obtained from a quintuple-chord occultation event involving Trans-Neptunian Object (174567) Varda (formally 2003 MW 12 ). The event was observed on September 10 th , 2018, and is the first-ever observed stellar occultation by this body; it is used to constrain Varda's size, 3D shape, geometric albedo, and bulk density. This campaign was carried out within the Lucky Star programme 1 , which aims at studying and characterising objects of the outer Solar System object, in particular Centaurs and Trans-Neptunian Objects (TNOs) in collaboration with the RECON project (Buie & Keller 2016).
Varda belongs to the hot classical TNO population. Moreover, it is a binary object, its satellite Ilmarë having about half of its size (Grundy et al. (2015). As such, the mass of the system can be derived from Kepler's law thus giving a direct access to Varda's bulk density, once the volume is pinned down using stellar occultations. Density measurements, in turn, contain clues of the conditions in the primordial solar system and/or the current internal structure of the body and other TNOs similar to Varda (Fernández 2020a; Barucci & Merlin 2020).
A review on Trans-Neptunian Binaries (TNB) is given by Noll et al. (2020). These authors outline the current understanding and progress in the study of TNBs from the identification of binaries, determination of their mutual orbits, measurements of the system's mass, density, rotational state, component colours, mutual events as well as the shapes. They report on 86 known TNBs at the time of writing. There are now 107 known binary and multiple systems in the TNO population (including Pluto and Charon) with one or more companions 2 .
This paper is organised as follows. Section 2 provides an overview of the Varda-Ilmarë system. Section 3 describes our occultation prediction approach. Section 4 presents the data, and the derived timings for each of the positive chords. In Section 5, Varda's projected size and geometric albedo are obtained and its 3D shape and density are constrained, and implications are discussed. Section 6 concludes the paper.

Varda's physical properties
Varda was discovered on June 21 st , 2003 with the 0.9-m Spacewatch telescope at Kitt Peak Observatory (Larsen et al. 2007). A satellite was discovered around it in 2009, see sub-Section 2.2. Varda's system physical properties, available in the literature, are summarised in Table 1. With an inclination of i = 21.5 • with respect to the Ecliptic plane, Varda falls into the category of the dynamically hot population, characterised by inclinations larger that 5 • (Fernández 2020b).
Following the new TNO taxonomy, which is based on colour indices (B-V, V-R, V-I, V-J, V-H, and V-K see Barucci et al. 2005;Fulchignoni et al. 2008), four classes are distinguished with increasingly red colours BB (neutral colour), BR, IR, and RR (very red). Using visible colour photometry, and applying this new taxonomy, Perna et al. (2010) have re-classified Varda as an IR object. After conducting a study on 75 known Centaurs and TNOs, Barucci et al. (2011) found that IR class objects belong to classical and resonant populations. None of the objects they have identified in this class, including Varda, seemed to contain an unambiguous water ice signatures near 1.5 µm and 2 µm. Furthermore, Varda showed the largest positive slope in the 2.05-2.3 µm range. The spectrum obtained by Barucci et al. (2011) is consistent with the presence of ice tholin on Varda's surface.

Varda's satellite, Ilmarë
The Hubble Space Telescope (HST) Snapshot programme 11113 (conducted between July 2007 and April 2009) allowed the observation of 142 TNOs with the Wide Field and Planetary Camera 2 (WFPC2). HST images obtained on April 26, 2009 revealed Varda's satellite, Ilmarë, at a separation of about 0.12 arcsec (corresponding to ∼ 4000 km, i.e. about 12 Varda radii), with a size of about half that of Varda (Table 1). Additional (visible and near IR) data, acquired between April 2009 and July 2013 with HST, Keck II Telescope, and Gemini North Telescope provided the relative offsets of Ilmarë with respect to Varda at twelve different epochs (Grundy et al. 2015 and references therein).
The derived Ilmarë's orbit has a mirror-ambiguity, two solutions being possible (cf. Table 1). The ambiguity should be removed by the end of 2020, when new relative positions will allow determining the correct solution, due to the changing aspect of the orbit as seen from Earth.
Furthermore, Grundy et al. (2015) used their HST observations to perform separate photometry of the Varda-Ilmarë system (translated in B, V, I Johnson magnitudes). The B-V magnitudes show that the components are consistent with one another, whereas the V-I magnitudes show an Ilmarë slightly redder than Varda (by (0.133 ± 0.062), see Grundy et al. 2015, and references therein). Moreover, Grundy et al. (2015) measured the difference in magnitude between the two components, ∆m = (1.734 ± 0.042) mag.  Vilenius et al. (2014) define the radiometric (area-equivalent) effective diameter of a binary system as D = D 2 1 + D 2 2 , where D 1 and D 2 are the primary's and secondary's diameters, respectively. (b) This assumes a linear trend of the phase curves and considers a magnitude variability that is due to the rotational light-curve. Thirouin et al. (2010Thirouin et al. ( , 2014 reported on broadband CCD photometric observations carried out between 2006 and 2013 using several telescopes for the Varda-Ilmarë system. They give three possible solutions for Varda's rotation period (Table 1). Although their analysis of the low-amplitude (0.06 ± 0.01 mag) single-peaked light-curve favours the 5.91 h period, the study shows the existence of two aliases at 4.76 h and 7.87 h that cannot be discarded. Note that all these periods are much smaller than the mutual orbit's period, which is estimated at ∼138 h (5.75 days), see Grundy et al. (2015). These authors compiled all the photometric data from Thirouin et al. (2010Thirouin et al. ( , 2014 and re-evaluated the associated Lomb periodogramme (Lomb 1976;Press et al. 1992), attempting to identify frequencies that could/would be expected for tidally locked components in the system. The results were inconclusive.

Varda's ephemeris & the occultation's prediction
Occultations by TNOs provide valuable opportunities to measure their physical characteristics from ground-based observations. Huge efforts are dedicated to the predictions of events (Assafin et al. 2012;Camargo et al. 2014;Desmars et al. 2015)). The second release of Gaia stellar catalogue (Gaia Collaboration et al. 2016 ensures an accuracy to the tenth of milli-arcsec (mas), thus most of the uncertainty in the prediction resides in the ephemerides of the TNOs.  Fig. 1: Difference between NIMAv6 and JPL7 in DEC (declination, left) and R.A. * (right ascension weighted by cos δ, right) over the time interval 2013-2022. Each dot represents the mean residual and the standard deviation of one night of observations. Dark blue points indicate data reduced with the WFI (Wide Field Imager) catalogue (see text), light blue points are for data reduced with Gaia-DR1 and light red points are for data reduced with Gaia-DR2. The continuous black line represents the best-fitting curve to these data, whereas the grey shaded area represents the 1-σ uncertainty of the NIMAv6 ephemeris (https://lesia.obspm.fr/luckystar/obj.php?p=46). To get the most accurate ephemeris, regular astrometric observations are performed by our team in different observatories (ESO, Observatório do Pico dos Dias, Sierra Nevada, etc). The astrometric reduction of these observations makes use of the Gaia catalogue, while orbits and ephemerides are obtained using the NIMA procedure (Desmars et al. 2015).
For the prediction of this occultation, astrometric observations from MPC (1980MPC ( -2018, ESO (2013), Observatório Pico dos Dias (2014Dias ( , 2017Dias ( , 2018 and Sierra Nevada (2018) were used to derive NIMAv6 solution of Varda's orbit. As the position of Ilmarë around Varda is unknown, photo-centre and barycentre of Varda's system are assumed to be blended with Varda's centre in the NIMA orbit determination. Figure 1 presents a comparison between the NIMAv6 and JPL7 solutions in right ascension (weighted by cos δ) and in declination between 2013 and 2022. Residuals of observations from our survey are represented by mean points with their standard deviations as error bars (one point for observations during the same night). Dark blue points indicate data reduced with the WFI catalogue (Assafin et al. 2012;Camargo et al. 2014), light blue points are for data reduced with Gaia-DR1 and light red points are for data reduced with Gaia-DR2. NIMA also allows an estimation of the orbit uncertainty (1-σ) represented in grey area.
Using NIMAv6 solution for the ephemeris and Gaia-DR2 stellar position, the predicted geocentric mid-time for the occultation was on September 10 th , 2018 at 03:38:37 UT (±117 sec). Table 2  This occultation was also predicted by RECON 3 using ephemeris based only on Minor Planet Center data without additional astrometric positions. The predicted path 4 , shifted to the West (∼ 1500 km) compared to Lucky Star prediction, highlights the importance of adding new astrometry and doing specific work to support the orbit fitting beyond automatic tools. Our new astrometric positions for Varda are provided in the appendix as well as a supplementary material in the form of an ASCII file.
The final prediction was made available to observers through the Lucky Star webpage, the RECON webpage as well as Occult Watcher 5 .  Table 2. Dark-grey, light-grey and clear areas are in the night, twilight and daytime, respectively. Varda's shadow is delimited by the continuous blue lines. It moves from north to south, each bullet corresponding to one-minute intervals, the largest bullet corresponding to geocentric closest approach. The red dashed lines show the 1-σ path uncertainty, largely dominated by the uncertainties of the NIMAv6 ephemeris ( Fig. 1) when compared to the star's position (Table 2).

Data analysis
Observations were performed at thirteen different stations in the USA (see acronyms in Table 3). The observers were organised in eight RECON teams as well as five IOTA 6 (non-RECON) teams. The BDR team used a QHY174M-GPS, while all the other stations used IOTA-VTI timestamp. All other relevant information regarding the observing circumstances are summarised in Table 3.

Occultation light-curves
Eleven stations recorded data during the occultation, five of which (TWF, MHV, YMA, CAR, and FLO) reported positive detections. Differential aperture photometry was used to obtain the occultation light-curves, using the Platform for Reduction of Astronomical Images Automatically (PRAIA, Assafin et al. 2011). The light-curves from the five aforementioned stations are shown in Fig. 3; their photometric measurements are provided as supplementary material to this paper. Fig. 3 presents these light-curves to which we have applied a manual offset in flux and centred at their respective mid-occultation times, for the sake of legibility.

Ingress and egress timings
The start (ingress) and end (egress) times of the occultation were determined by fitting each event by a sharp opaque edge model, after convolution by (1) Fresnel diffraction, (2) finite CCD bandwidth, (3) finite stellar diameter, and (4) finite integration time (see Braga-Ribas et al. 2013;Ortiz et al. 2017 for details). 3 Research and Education Collaborative Occultation Network (Buie & Keller 2016) 4 https://www.boulder.swri.edu/~buie/recon/events/174567_180910_0154895.html 5 https://www.occultwatcher.net/ 6 The International Occultation Timing Association has had a four decade old mission to observe and record eclipses of stars by minor planets and other bodies. Since the first announced observation there have been thousands of successful observations. In the past decade digital video records have been created with time sensitivities to .003 seconds. IOTA contributors did not receive funding and participated in this study on a voluntary basis.
(a) Post-occultation map.  The red pins/ black pins on the map represent the sites that reported +ve / -ve observations, whereas the white ones are for the (NDR -no data recorded) stations that experienced technical issues preventing them from recording the occultation (cf. Table 3). Sub- figure 3b shows the five occultations detected at the TWF, MHV, YMA, CAR, and FLO stations (from the western-most to the eastern-most on the Sky-plane cf. Figure 6). Each light-curve has been normalised between 0 (Varda's flux) and 1 (unocculted star + Varda's fluxes). For the sake of legibility, each light-curve has been centred at mid-chord time for each chord and a manual offset in flux has been applied. The associated data to the normalised, non-centred light curves is provided as supplementary material to the paper.
Fresnel diffraction operates over the scale F = √ λ∆/2, where λ is the wavelength of observation and ∆ the geocentric distance of the object. From Table 2 we obtain F ∼ 0.5 km in visible wavelengths. We estimate the stellar diameter to 0.8 km projected at Varda's distance, using B, V, K magnitudes of 15.61, 14.32, 11.704 for the star (NOMAD catalogue, Zacharias et al. (2004)), respectively, using van Belle (1999)'s formulae and Ochsenbein et al. (2000) 7 . Given the geocentric shadow velocity of 10.86 km s −1 and exposure times > 1 s at all stations (see Tables 2 and 3), the main cause of light-curve smoothing is the finite integration time, and not Fresnel diffraction or stellar diameter. However, the CAR station required further analysis, as an instrumental effect caused a smoothing of the light-curves over several data points (see Fig 5), resulting in gradual ingress and egress, see discussion below. The occultation timings are obtained by minimising the classical χ 2 function, using the same procedures as described in Sicardy et al. (2011). The best-fitting functions are displayed in Fig. 4 and the resulting chord details are listed in Table 4.

The case of the CAR station
The CAR station light-curve (Fig. 5) shows a gradual ingress and egress of the star; both the ingress and egress extend over more than 5 data points, i.e. more than five seconds. Two scenarii were considered to explain this behaviour, (i) a Varda atmosphere or (ii) an instrumental effect.
In case (i), a ray-tracing code was considered. It uses a nitrogen atmosphere with a temperature profile comparable to that of Pluto, with a rapid increase of temperature just above the surface that connects to an upper isothermal branch at ∼100 K. Then a surface pressure and temperature of ∼ 1 µbar and 38 K, respectively reproduces the profile of Fig. 5. Meanwhile, the four other stations (TWF, MHV, YMA, and FLO) do not have sufficient SNR and/or time resolution to test this model. Case (ii): a closer examination of the CAR images revealed a slow motion of the stellar images. During ingress, the occulted star fades away but its image no longer moves on the array, while the other sources continue to move clearly, indicating an instrumental effect (more details are given in Appendix B). To account for the aforementioned effect, we have added in the modelling of the occultation light-curve a non-instantaneous instrumental response. More precisely, we have considered two possibilities for the response to an instantaneous light pulse: a. in the first case, we assume a simplified response that decays linearly over a time interval of ∆t resp . Exploring values of ∆t resp between 0 and 10 seconds with steps of 0.01s, the occultation timings and the optimal ∆t resp values timings are obtained by minimising the classical χ 2 function as done for the other light-curves. We find the best (and satisfactory) simultaneous fit to the ingress and egress for ∆t resp = 6.58 +0.63 −0.52 s. (Fig. 5a).  b. in the second case, we assume an (RC) electronic filter response (e −t /τ , where τ is the capacitor's time constant) convolved with a rectangular function response that accounts for the finite integration time (cf. Table 3). We explore τ values in the interval [0,10] seconds with a step of 0.015 s. The occultation timings and the optimal τ values timings are obtained by minimising the classical χ 2 function as done for the other light-curves. We find the best (and satisfactory) simultaneous fit to both ingress and egress for τ = 3.21 ± 0.45 s (Fig 5b).
The resulting ingress and egress times obtained by the two methods agree at the 1-σ confidence level. However, further analysis of the post-occultation data (cf. Appendix B) motivated our choice of an electrical RC filter (case b above). At last, the overall uncertainty in the timing of this chord is mainly (∼ 60%) due to the uncertainty in the modelling of the (unknown) instrumental response. The gradual ingress and egress is clearly due to an instrumental effect.
The resulting egress and ingress timings as well as the duration and length of the chords at the five positive stations are listed in Table 4.

Elliptical fit to Varda's limb
Given that Varda's diameter is larger than 500 km (cf. Table 1), it is reasonable to assume that its limb is close to elliptical due to hydrostatic equilibrium (Tancredi & Favre 2008). Therefore, following the same procedures as in previous works (   i. The centre of the ellipse, ( f c , g c ), that measures the offsets in α and δ to apply (eastwards and northwards, respectively) to the adopted ephemeris (NIMAv6), assuming the star position of Table 2; ii. The apparent semi-major axis a ; iii. The apparent oblateness = (a − c )/a , where c is the apparent semi-minor axis; iv. The position angle P of the apparent semi-minor axis c , measured eastwards from Celestial North.
These parameters are varied in wide ranges to find the best-fitting limb. More precisely, we explore: −90 • < P < 90 • with steps of 1 • , 0 < < 0.5 with steps of 0.005, and 350 km < a < 450 km with 1 km steps. For each value of P , and a , we determine the best-fitting ellipse by adjusting the centre ( f c , g c ). This best-fitting limb is obtained by minimising the function χ 2 = N=6 i=1 (r i − r ell,i ) 2 /σ 2 i , where r i , r ell,i the radial distance of the i th occultation point to the shadow's centre and to the elliptical limb model, respectively. σ i is the radial uncertainty (1σ level) stemming from the occultation's timing uncertainties (Table 4). The χ 2 value is based on the radial residual of the actual limb relative to the limb model. The associated error bar on each data point is then the timing uncertainty in the direction of the chord multiplied by the (radial) velocity of the star with respect to the centre of the body. As such, it is fully consistent with the approach that would minimise the residuals along the chords. Admittedly, this equivalence of methods works as long as the limb remains close to circular and none of the chords are close to grazing, which is the case here. We define the value of χ 2 per degree of freedom (or unbiased χ 2 ) by χ 2 pdf = χ 2 /(N − M), where N = 10 is the number of data points and M = 5 is the number of adjustable parameters. Fig. 6 displays the results of those fits and Table 5 provides the best-fitting adjusted parameters. The expected minimum value of χ 2 pdf for a satisfactory fit is χ 2 min pdf = N−M M = 1, which is the case from Table 5, noted χ 2 pdf = 1.32.  Table 4, while the red segments are the uncertainties at the chords extremities. The black ellipse is the best-fitting solution, whose parameters are given in Table 5. The grey area shows all the solutions that are found to within a 1-σ level from the best fitting ellipse, i.e. corresponding to χ 2 < χ 2 min + 1. The purple dashed lines show the projections onto the sky from the two NDR stations (cf. Table 3), whereas the brown dotted lines show the projections from the near miss stations (cf. Figure 3a). Notes. (a) The apparent area-equivalent radius is R equiv = √ a c = a √ 1 − , where c is the apparent polar radius. (b) We do not know Varda's rotational phase at the epoch of the occultation; therefore we cannot account for the peak-to-peak ∆mag = 0.06 mag amplitude of Varda's rotational light-curve (Thirouin et al. 2010). This adds an additional uncertainty of ±0.03 on p. (c) The associated centre ( f c , g c )=(50.49±3.29, 59.43±1.81) km of the ellipse (Fig. 6) measures the offsets in α and δ to apply to the adopted NIMAv6 ephemeris, which yields to this new astrometric position. (d) The calculation of χ 2 pdf is explained in the text.
In our paper, the 1-σ error bars on each fitted parameter stem from their respective marginal probabilities, and more precisely, their 68.3% confidence level values, without consideration of the other adjusted parameters. These domains are bound by the χ 2 < χ 2 min + 1 criterion. This same approach was adopted in (Braga-Ribas et al. 2013Benedetti-Rossi et al. 2016 Article number, page 9 of 17 Dias-Oliveira et al. 2017;Ortiz et al. 2017). The χ 2 < χ 2 min + 1 limit corresponds to confidence limit for one free parameter, while considering simultaneously all five free parameters, the formal 1 − σ confidence value for χ 2 would write χ 2 < χ 2 min + 5.89 (this is not the adopted approach here).

The geometric albedo p v
We use the equivalent radius R equiv given in Table 5 and Varda's visual absolute magnitude H V from the literature to derive the associated geometric albedo for Varda. The absolute magnitude given in the literature is that of the system Varda-Ilmarë (cf . Table 1). Moreover, as mentioned in Sect. 2.2, Grundy et al. (2015) have measured a difference in the visible magnitude between the two components ∆m = (1.734 ± 0.042) mag. It follows that Varda's absolute magnitude is H v = 3.81 ± 0.01. Using this latter value of the absolute magnitude, we derive the geometric albedo of Varda 8 , p = 0.097 ± 0.004.
The geometric albedo value that we derive is consistent with but more accurate than the value from radiometric measurements p V = 0.102 +0.024 −0.020 of Vilenius et al. (2014). This shows that Varda is about as dark as the binary TNO 2003 AZ 84 (Dias-Oliveira et al. 2017).
Moreover, recent thermal emission data of the Varda-Ilmarë system observed with ALMA has shown that the ratio of the primary's radius to the secondary's is ∼ 2.10 (Moullet, A. & Lellouch, E., priv. communication May 2020). Given the equivalent radius of Varda (cf. Table 5), we derive an equivalent radius of ∼ 178 km for Ilmarë, about 15 km larger than Grundy et al. (2015)'s value. Furthermore, using a difference ∆m = 1.734 mag between the primary and the secondary implies that the albedo of Ilmarë is ∼0.89 times that of Varda, i.e. ∼0.085.

Geometrical considerations
We now use the 2D limb fit obtained above to constrain Varda's 3D shape. To do so, simplified assumptions must be made: Varda is an oblate, axisymmetric spheroid with equatorial radius (a), polar radius (c), and true oblateness = (a − c)/a. This choice is based on the small amplitude of Varda's rotational light-curve (∆mag = (0.06 ± 0.01) mag. 9 , see Table 1) which favours a spheroid solution over a non-axisymmetric, triaxial ellipsoid (Thirouin et al. 2014), i.e. the flux variations caused by albedo variegation.
The axisymmetric assumption leads to a = a . Moreover, the apparent oblateness is related to the true oblateness by the classical formula c 2 = a 2 sin 2 B + c 2 cos 2 B (c is the apparent polar radius), i.e.
where B is the planeto-centric sub-observer latitude, such that B = 0 • (B = 90 • , resp.) corresponds to equator-on (pole-on, resp.) viewing geometry. The density ρ is then derived from the mass M divided by the volume 3a 2 c/4π. From c = a(1 − ) and Eq. 1, we obtain The mass M is an observable quantity, while a and are derived in the previous sub-section (see Table 5). Grundy et al. 2015 give the system's mass (Varda plus Ilmarë) M = (2.664 ± 0.064) × 10 22 kg, where Varda accounts for 92% of the system's mass. We adopt here a mass of M = (2.45±0.06)10 22 kg for Varda alone. As B is a priori unknown, Eqs. 1 and 2 define a parametric curve (ρ) travelled as B varies from 0 (the equator-on geometry with minimum possible value of = ), to its maximum value arcsin(1 − ), corresponding to a completely (and unrealistic) flattened object, i.e. = 1. The function (ρ) is plotted in Fig. 7 as the light-blue dashed line. The uncertainty in causes a displacement of a given point ρ( ) along this same curve, while the uncertainties in M and a cause a transverse relative uncertainty in ρ of: since M and a are determined independently. The resulting error domain for ρ is delimited by the light-blue continuous lines in Fig. 7. The minimum density provided by the spheroid model is ρ ∼ 1.16±0.04 g cm −3 , up to an infinite value at B = arcsin (1− ) ∼ 73 • , corresponding to a flat object with = 1.

MacLaurin solutions
We now go a step further, assuming that Varda is homogeneous and in hydrostatic equilibrium. This implies that Varda is either a MacLaurin spheroid, or a Jacobi triaxial ellipsoid (Chandrasekhar 1987). Note that with a typical equatorial radius of more than 380 km (Table 5) and from its derived density (see below), Varda qualifies as a candidate for being a dwarf planet, i.e. a body in hydrostatic equilibrium (Tancredi & Favre 2008). As mentioned earlier, an axisymmetric spheroid (here, a MacLaurin solution) is preferred over a triaxial ellipsoid (a Jacobi solution). In this case, the density ρ of the body is related to its spin frequency ω = 2π/T rot (where T rot is the rotational period) and its true oblateness by where cos(ψ) = 1 − and G is the gravitational constant (Plummer 1919;Chandrasekhar 1987;Sicardy et al. 2011). Three possible rotational periods (for single-peaked light-curves) are reported in the literature (see Section 2.3): 4.76 h, 5.91 h (the most probable, according to the authors), and 7.87 h. The resulting MacLaurin curves vs. ρ are shown in Fig. 7 as purple, green and red curves, respectively. The intersections of the light-blue curve and the MacLaurin solutions in Fig. 7 provide solutions that are consistent with both our occultation results and the MacLaurin hypothesis. All possible MacLaurin solutions consistent with our occultation's results are summarised in Table 6. Solutions 1, 2, and 3 are associated with the 4.71 h, 5.91 h, and 7.87 h rotational periods of the single-peaked light-curves, respectively; while solutions 1bis, 2bis, and 3bis are associated with the respective double periods (i.e. 9.52 h, 11.81 h, and 15.74 h). Although the light-curves associated with these double periods would mean a body with a double-peaked rotational light-curve due to albedo features on opposite sides of the body, these solutions shall not be discarded because of the low-amplitude of the rotational light-curve (see previous sections). This stems from the extremely low-amplitude variability (∆m < 0.15 mag), for which distinguishing the period from its double is nearly impossible (Thirouin et al. 2014). Finally, we cannot exclude the existence of a non-axisymmetric, triaxial ellipsoid solution, in particular a Jacobi-type body close to a pole-on orientation (to explain the low peak-to-peak amplitude of the light-curve ∼0.06 mag). Asides for solution 1 (for which T = 4.76 h), all the others give bulk density values that are consistent with but more accurate than the value 1.27 +0.41 −0.44 g cm −3 derived from thermal measurements for which equal albedo values and equal densities were assumed for both Varda and Ilmarë (Thirouin et al. 2014). Solution 1 gives a much higher and seemingly unlikely density for this size of object (2.18 ± 0.07 g cm −3 ), close that of Quaoar. This solution, however, shall not be excluded (see discussion in Sect. 5.3).
A further constraint can be considered by requiring an equatorial orbit for Ilmarë. This may, for example, be expected if Ilmarë resulted from the re-accretion of a collisional disc surrounding Varda after an impact. The opening angle of its orbit, B open , can be calculated from the orbit's pole position (Table 1) and Varda's astrometric position at epoch.
We examine here the possibility that Ilmarë moves in the equatorial plane of one of the MacLaurin solutions considered above. Grundy et al. (2015) provide two (mirror-ambiguous) solutions for Ilmarë's orbital pole (Table 1) Fig. 7. Moreover, it would require an unrealistic high density for Varda, ρ > 2.5 g cm −3 . In contrast, the second Ilmarë orbital solution (represented by the shaded area in Fig. 7) lies between solutions 2 and 3 in terms of opening angles. Moreover, it provides a position angle for the semi-minor axis of the orbit (as projected in the sky plane) of 96 • .2 +4.0 −3.6 . This is consistent with the position angle P = (113 • ± 34 • ) of Varda's limb semi-axis (Table 5), and thus, with an equatorial orbit for Ilmarë. This solution has a bulk density and oblateness in the ranges 1.47-1.71 g cm −3 and 0.16-0.20, respectively (see the light-blue shaded region in Fig. 7).

Discussion
The geometric albedo we derive for Varda (p = 0.097±0.004, see Table 5) can be compared to those of other TNOs. The dynamically hot population has a median geometric albedo of 0.085 +0.084 −0.045 (Müller et al. 2020), excluding the Haumea family 10 and dwarf planets, noting that the median geometric albedo for cold classical is 0.14 +0.09 −0.07 (Ibid.). More recently, Müller et al. (2020) provided an overview of TNO physical properties at thermal wavelengths, in particular, for a sample of 26 hot classical TNOs (excluding Makemake, Haumea and its family), with geometric albedo values ranging between 0.032 and 0.310, a median at 0.084, and a mean value at 0.102. Thus, Varda's geometric albedo appears in line with those of other hot classical TNOs.
Our results are based on the assumptions that Varda is an oblate, axisymmetric spheroid with equatorial radius (a), polar radius (c), and true oblateness = (a−c)/a. The small amplitude of Varda's rotational single-peaked light-curve (∆mag = 0.06±0.01 mag., see Table 1) favours a spheroid solution over a non-axisymmetric, triaxial ellipsoid, and motivates our search for MacLaurin solutions.
With an equatorial radius of about 380 km, Varda falls into the category of Kuiper Belt Objects (KBOs) that make a transition from small porous objects to dense KBOs, as defined in Bierson & Nimmo (2019). Simulations by those authors show that such a transition rapidly occurs for objects in the radius range 200-500 km. To explain the observed density distribution of the KBOs as a function of size, Bierson & Nimmo (2019) studied a sample of 11 out the 18 well-known TNBs using a 1D model that couples KBO's thermal evolution with porosity evolution (where primordial porosity is removed over time). They have concluded that the density distributions observed within the KBO populations are mainly a consequence of porosity, rather than mass. They have also defined an input rock mass fraction parameter 11 f m = M S /(M S + M i ), where M S is the mass of silicates and M i the mass of ice. In their model Varda has f m > 70%, close to Pluto's value, and much lower than those of Eris ( f m ∼ 90%) and Quaoar ( f m ∼ 85%) which are much larger and much denser objects compared to Varda. This confirms that the Solution 1 given (ρ = 2.18g cm −3 ) in Table 6 cannot be excluded just yet.
To summarise, solutions 1, 2, and 3 are compatible with a spheroid MacLaurin solution, whereas solution 1bis, 2bis, and 3bis would be associated with triaxial ellipsoid shape of Varda. As discussed in sub-section 5.2.3, distinguishing a period from its double is nearly impossible when dealing with low-amplitude light-curves. As an example, one can cite Ceres, whose oblate shape is wellknown from stellar occultations (Gomes-Júnior et al. 2015) and the DAWN spacecraft visit (Russell et al. 2016), while it exhibits a low-amplitude double-peaked rotational light curve due to albedo features (Chamberlain et al. 2007). In the trans-neptunian region, Makemake is also a clear example of oblate body with low-amplitude double-peaked light-curve (Hromakina et al. 2019).
If we were to assume Varda as a body with a double-peaked light-curve due to albedo features on opposite sides of the body, we should also consider the double of the aforementioned periods (i.e. 9.52 h, 11.82 h, and 17.74 h, respectively). The associated MacLaurin curves are plotted in figure 7 in purple, blue, and red dashed lines, respectively. The associated density and oblateness values are given in Table 6.

Conclusion
This study provides unique constraints on fundamental physical properties of Varda, by refining its size, shape, geometric albedo, and bulk density. It uses the first-ever recorded stellar occultation by the TNO (174567) Varda, observed from several stations in the US on September 10 th , 2018. This event resulted in five positive chords, from which an elliptical limb-fitting provides an apparent semi-major axis of (381 ± 3) km and an apparent oblateness of 0.043 ± 0.036, respectively, an area-equivalent radius of 373 ± 8 km, and a geometric albedo of p v = 0.097 ± 0.004.
We have derived six possible MacLaurin solutions for Varda, three assuming single-peaked rotational light curves and three assuming the associated respective double-peaked rotational light-curves. The associated density and oblateness values are listed in Table 6.
At this time, we cannot discriminate between these solutions. However, two solutions require our attention: that associated with the most probable period (T=5.91 h) for a spheroid body with a single-peaked light-curve (solution 2: ρ = 1.52 ± 0.05 g cm −3 and = 0.232 ± 0.036), and its double period (T=11.82 h) for a body with a double-peaked light-curve (solution 2bis: ρ = 1.25 ± 0.04 g cm −3 and = 0.079 ± 0.044). The relatively high densities derived for Varda (>1.5 g cm −3 ) adopting the McLaurin solutions 1, 2, and 3 are consistent with low porosity for an object in this size-range.
The gradual ingress and egress observed at one of the stations mimic the effect of an atmosphere, but are actually of an instrumental origin. Such effects could be observed again when using similar setups, so some care is necessary to avoid erroneous interpretations.
From a more general standpoint, Trans-Neptunian Binaries like Varda are good occultation candidates to probe the outer Solar System, because they provide accurate shapes (i.e. volumes) for bodies whose masses are well constrained from the satellite's motion. They provide a wide range of densities, from below that of water ice to that of nearly pure rock. Such occultations constrain the formation scenarii for the Solar System and/or the current internal structure of those bodies.

Acknowledgment
This campaign was carried out within the "Lucky Star" umbrella that agglomerates the efforts of the Paris, Granada and Rio teams. It is funded by the European Research Council under the European Community's H2020 (2014-2020/ERC Grant Agreement No. 669416) . Fig. B.1: This figure shows a sequence of 8 sub-frame images taken a few minutes after the event (between 03:43:54 and 03:44:02 UT). On each frame, the target star is circled in red while the guide and calibration stars are circled in yellow and green, respectively. Due to the tracking correction (because of imprecise telescope tracking), the stars move abruptly on the CCD during the exposure at 03:44:55 (on the sub-frame, we only see the minutes and the seconds). The green arrows point towards the real stars' positions. On the 03:43:57 exposure for example, we see "phantoms" of the dimmest stars (remnants from the previous exposures), which are shown by the white arrows. These phantoms are seen for at least five images, up to the 03:44:01 exposure.
Article number, page 15 of 17 A&A proofs: manuscript no. Souami_et_al  Fig. B.3: This figure shows 8 successive integrations of 1.28 s with a Watec 910HX and 3DNR set to ON at 50% level, during an abrupt correction of the star field centring. The green arrows point towards several faint reference stars. The white arrows point towards ghost images of these faint stars at their first position, fading gradually with time. Bright stars are not impacted by the filter.