Stellar Occultation by the Resonant Trans-Neptunian Object (523764) 2014 WC510 Reveals a Close Binary TNO

We report on the stellar occultation by (523764) 2014 WC510 observed on 2018 December 1 UT. This occultation campaign was part of the Research and Education Collaborative Occultation Network (RECON), a network of small telescopes spread over 2000 km in western USA and Canada. Light curves from six stations revealed three groups of two or more consecutive flux drops correlated in time between adjacent stations. A Bayesian model comparison reveals that a model with a double object occulting a double star is favored over alternative models considered. For the statistically favored model, we determined that the primary component of the object has a diameter dp = 181 ± 16 km and the secondary ds = 138 ± 32 km, assuming identical geometric albedo between the two components. The two components have a projected separation of 349 ± 26 km. Adopting an absolute magnitude for the system of HV = 7.2 from the Minor Planet Center, we derive a geometric albedo of pV = 5.1% ± 1.7%. This is the smallest resonant object with an occultation size measurement and with a detected secondary from a ground-based stellar occultation, filling a region of the size versus separation parameter space of binary objects that is largely unexplored. The results show the capabilities of the unique design of the RECON experiment sensitive to small objects and close binaries. 2014 WC510 is presently at a low galactic latitude where the high surface density of stars will provide good occultation opportunities in the upcoming years.


Introduction
(523764) 2014 WC 510 was discovered on 2011 by Pan-STARRS1 (Chambers et al. 2016). The Deep Ecliptic Survey (DES) dynamical classification of this object is that it inhabits a mean-motion resonance with Neptune, specifically the 3:2e + 6:3i 2 . See Elliot et al. (2005) for details about the DES classification system. This object is also in a Kozai secular resonance with an 18°libration amplitude. In comparison, the Kozai libration amplitude of (134340) Pluto is 24°.4. However, the Kozai resonant argument for 2014 WC 510 is not as cleanly sinusoidal as that for Pluto. The orbital period is 246 yr, its heliocentric distance ranges from 28.7 to 50.8 au, and its orbital inclination is 19°.8.
Very little is known so far about this object's physical properties. At the time of the occultation, it had an apparent magnitude of V=22.1 and was 30.5 au from the Sun and 29.6 au from Earth. Based on its estimated absolute magnitude of H V =7.2±0.3, 2014 WC 510 would have a diameter of 220 km assuming a 5% albedo, or a diameter of 90 km assuming a 30% albedo.
Size, geometric albedo, and binarity are basic physical properties of trans-Neptunian objects (TNOs), and these are the primary observational goals of the Research and Education Collaborative Occultation Network (RECON). With these measurements, we expect to contribute to the knowledge of the size distribution of these objects. Additionally, the binary fraction among different dynamical classes and binary properties such as component size ratio and inclination give us clues about the prevalent planetesimal formation mechanism . Furthermore, binaries with known orbits can provide direct measurement of the system mass and ultimately the bulk density when combined with size measurements. Stellar occultations provide a powerful tool to support the collection of these measurements.
An estimate for the size of TNOs can also be obtained from radiometric measurements with typical uncertainties of 10%-20% (Moullet et al. 2011). More than a hundred objects have size estimations from this technique (Kovalenko et al. 2017). The stellar occultation technique provides more accurate measurements of sizes and shapes, but the thermal method, so far, is far more efficient in the collection of large numbers of measurements, albeit of lower precision. Sizes from stellar occultations are directly useful but can also provide crucial data for cross-calibration of thermal measurements . To date, only 14 TNOs or their satellites (excluding Pluto) have an accurate size and albedo published from stellar occultations: (134340) Pluto I Charon (satellite of Pluto; Sicardy et al. 2006 (Ortiz et al. 2020), (541132) Leleākūhonua (Buie et al. 2020a), and (486958) Arrokoth (Buie et al. 2020b). Those with successful occultations are mainly limited to the largest objects (D200 km) because those are easier to predict and thus observe. With the availability of accurate star astrometric positions provided by the Gaia mission, currently in its second data release (Gaia DR2), it is now feasible to attempt stellar occultations of smaller objects (Porter et al. 2018;Buie et al. 2020b). Even with improved predictions, successful occultation results demand a larger effort than is typically needed with larger objects. The stellar occultation by the TNO (229762) G! kún'hòmdímà represents an excellent case of a complete characterization of a binary TNO, measured with an occultation in 2014 by RECON (Benedetti-Rossi et al. 2016), and its orbit later characterized by astrometry from the Hubble Space Telescope (HST; Grundy et al. 2019b) .
Detection and characterization of binaries have largely been conducted with direct imaging with ground-based adaptiveoptics-assisted facilities and space-based observations from HST (Grundy et al. 2019a). The most notable exception is the contact binary cold-classical Arrokoth studied remotely by occultation (Buie et al. 2020b) and in situ by the New Horizons extended mission (Stern et al. 2019). Setting aside Arrokoth, the properties of the binary population are limited by the spatial resolution of direct images and as a result are limited to separations larger than ∼1000 km. The Arrokoth occultation result is a powerful example of how occultations probe inside this limitation of direct imaging.
Recent works based on rotational light curves hint that the number of contact binaries among the 3:2 resonant population is larger than those in the classical Kuiper-Edgeworth belt (Thirouin & Sheppard 2018). Only through stellar occultations can we hope to provide firm confirmation of these candidates and determine the binary or contact binary properties. RECON (Buie & Keller 2016) was designed for this purpose and can detect objects down to about 50-100 km in diameter while also searching for multiple systems with separations of about 2000 km all the way down to contact binaries.
Most stellar occultation results published so far involve large objects with relatively good signal-to-noise ratio (S/N) where the standard analysis approach is to extract the times of disappearance and reappearance of the star as seen from each observer. The determination of a physical parameter such as size, shape, and pole position is then determined by fitting a 2D or 3D shape using a minimization approach (Widemann et al. 2009) or using a Bayesian approach (Brown 2013;Leiva et al. 2017). The determination of disappearance and reappearance times becomes impractical when dealing with stellar occultations where the disappearance and reappearance of the occulted star are ambiguous because of low S/N, short duration of the occultation, or both. Such cases often require comparison of different models, including the case where the stellar flux might be solely explained by random noise variations in the light curve with no occultation at all.
In our work, we perform two levels of Bayesian inference (1) to quantitatively compare alternative models to determine which one is favored given the data and a model for the uncertainties in the measurements, and (2) to estimate the parameter values and uncertainties for the favored model. In the sections that follow, we provide information about the prediction of the occultation, the general conditions during the occultation campaign, the observations, and the data reduction. We will further provide the details of the data analysis, the different models considered and why, and our approach to comparing different models while deriving the nominal parameter values for the favored model and their uncertainties.

Occultation Prediction and Circumstances
The prediction and observations for this occultation campaign were part of the RECON project (Buie & Keller 2016). The prediction was done in the same way as for all other objects targeted by RECON. For each object, the RECON automatic prediction system uses all the astrometric individual positions available from the Minor Planet Center (MPC) to update the orbit of the object and to keep a daily updated list of favorable occultation events visible anywhere on Earth with a targeted list extracted that is specifically observable from the network. An occultation opportunity by 2014 WC 510 was automatically identified by the RECON prediction system with a geocentric close approach at 2018 December 1 12:21:23 UTC.
The occultation by 2014 WC 510 presents an interesting case. At the time of the event, all of the astrometric positions available at the MPC for this object were submitted by Pan-STARRS. The RECON project worked out a collaboration with the Pan-STARRS project through R. Weryk to mine the Pan-STARRS image archive for all the TNO observations not found through regular processing for near-Earth objects. Any such observations found are noted, measured, and reported to the MPC, which then soon flow into the RECON prediction system. This effort also included an upgrade in the Pan-STARRS astrometric process from using the 2MASS catalog to using the Gaia DR2 catalog (Brown et al. 2018). After this process, about 80% of the TNOs' astrometric positions submitted to the MPC by Pan-STARRS were reduced against Gaia DR2. By 2018 September 7, the number of observations for 2014 WC 510 was up to 89 and spanned a total arc length of 6.3 yr. Normally, an arc length this short is not good enough for a useful occultation prediction, but having data from a single source with an exceptionally good ground-based astrometric system was extremely effective. Table 1 summarizes the positional parameters of the occulted star obtained from the Gaia DR2 catalog. The astrometric position used in the prediction and occultation analysis is obtained from Gaia DR2 values using a linear approximation for the star's space motion and an approximate correction for the annual parallax with the resulting values and propagated uncertainties given in the lower section of Table 1. This includes the estimated systematic uncertainties for the parallax and proper motion from Lindegren et al. (2018), which we add in quadrature.
Details of the occultation prediction and observation circumstances are summarized in Table 2. The absolute magnitude in the V band is obtained from the MPC, while the apparent magnitude in the V band and the object-Earth distance are obtained from the RECON ephemerides system, which uses all MPC astrometric individual measurements (see Buie & Keller 2016 for details). From the same ephemerides, the sky-plane scale and geocentric occultation shadow speed at the predicted closest approach time are given, which must be taken only as informative values. The shadow speed changes as a function of time from the nominal value in Table 2 by a few meters per second, and it changes as well for each site. For the analysis, the informed values are not used, and instead the midtime of each data point in the light curves is transformed into a position in the sky plane. The sky plane is a plane perpendicular to the Earth-star line at the distance of the object with coordinate axes (x, y) in the direction of the east and north and with its origin coincident with the ephemeris of the center of the object (Elliot et al. 1978).
The target star was easily observable in dark skies for the entire RECON coverage area, with most sites seeing the star at an elevation between 30°and 40°. The cross-track uncertainty was 1100 km, and the event time uncertainty was 86 s (all uncertainties in this work are 1σ unless otherwise noted). The cross-track uncertainty in this case was dominated by the uncertainty in the orbit of the object with a comparatively negligible contribution from the uncertainty in the star position. Based on the nominal cross-track positions of the stations and the shadow-track uncertainty, there was a 35% chance of success for this event assuming all stations participated under clear skies.

Observations
The RECON stations involved in the observations are summarized in Table 3. From the 41 participating stations, 12 provided useful data with 6 of those (Bishop, CPSLO, Searchlight, Mohave Valley, Wildwood, and Parker) showing two or more consecutive flux drops below 2σ level from the light-curve average value. From the remaining sites, 26 had bad weather, one had technical issues during data transfer and the data were lost, one was unable to observe due to unavoidable obstruction of the field by the observatory structure and one had telescope alignment problems leading to observing the wrong field. In Table 3, the siteID is listed as described in Buie & Keller (2016). Sites with "C" codes correspond to the RECON 100 km Canadian extension CanCON incorporated in mid-2018. If data were collected, the start and stop times of collection are listed. For cases with no data, the nominal team location is listed, but for those sites with data, the actual GPSbased observing location is tabulated. Times are given in UTC while the locations are given in the WGS84 reference system and datum. The S/N column indicates the average signal-tonoise ratio of the light curve, determined as the ratio between the average signal and the standard deviation (see Section 4 for details in data reduction). The observers involved are listed followed by relevant comments regarding the data and weather conditions reported by the observers. For teams that were    Table 3. Superimposed on the map is the nominal predicted ground track, shown with a gray shaded and hatched region nearly centered on the Yuma site that is based on a diameter of 220 km corresponding to 5% geometric albedo.
The northern 1σ limit of the prediction is shown by the black dashed line passing between Bend and Chiloquin. The other elements in this figure will be addressed later in Section 6.  Table 3 and Figure 3 for more information about the sites and their data). The gray hatched region shows the RECON prediction using a diameter of 220 km (5% albedo) with the 1σ uncertainty in black dashed line. The red shaded regions are the shadow of the primary component (bottom) and secondary component (top) of the object occulting the primary star using the nominal parameter values of the favored model ( Table 7).

Data Reduction
All the data collected came from standard RECON sites, each consisting of an 11 inch telescope equipped with an interline-transfer CCD MallinCam camera used without filters, giving a field of view of 17′ by 13′ (1 6 px −1 ) and negligible dead time between images (∼1 μs). The raw data consists of video files with the GPS-based UTC time overlaid in the video; see Buie & Keller (2016) for a detailed description of the RECON video capture system.
The video files were converted to individual image files after performing a pixel-by-pixel robust average of all the video frames that belong to a single integration. The time of the start of each integration was derived from the time overlay inserted in the video frames as detailed in Buie & Keller (2016). Images were calibrated with dark and sky flats produced from short dark and flat calibration videos of about one-minute duration taken immediately after the occultation data.
The target star photometry was ratioed to the measurements of nearby, similar-brightness stars to remove trends produced by variations in transparency. This step was also useful to help discriminate false flux drops in the target star produced by image artifacts. Saturated hot pixels were the most common cause of such flux drops. The light-curve flux of the occulted star plus occulting object were then normalized to unity when the star is not occulted. In practice, the apparent magnitude of the occulting object (V=22.1) is much fainter than the occulting star (G=15.4), making the object flux negligible with respect to the flux of the star. At the same time, the apparent magnitude of the object is much fainter than the limiting magnitude of the images (V∼17), and the flux of the object is below the noise level of the images. The results for all sites with usable data are shown in Figure 3, highlighting in red all those data points at least 2σ away from the average lightcurve value. Error bars are omitted in the figure for clarity. The data used to construct this figure, including error bars, are provided as machine-readable tables.
The data from the (V-07) Wildwood site required some extra handling. The data were affected by abnormal camera operation. The team reported using a camera setup SEN-SEUP=128 (equivalent to an exposure time of 2.14 s) for data acquisition, but analysis of the data was consistent with SENSEUP=110 instead (equivalent to an exposure time of 1.83 s). This apparent integration time is not one of the valid settings for the camera software although the detector, being a CCD, is perfectly capable of this integration time. Missing frames would be evident in the video stream generated by the camera as jumps in the times saved in the video overlays (see Buie & Keller 2016 for a detailed description of the video capture system). We checked manually the time stamp overlay in each video frame to corroborate that there were no missing frames and that 1.83 s was the actual integration time of the camera. Additionally, this video suffered from an excessive noise pattern affecting about 45 rows in each video frame (about 10% of the frame). The affected rows moved vertically from frame to frame, thus affecting a different section in consecutive frames, allowing us to retrieve a cleaner image for the integration after performing the pixel-by-pixel robust average.

Data Analysis
The data were acquired with relatively long exposure times, either 2.14 or 1.83 s, corresponding to 48 km and 41 km in the sky plane, respectively. Given the Earth-object distance at the moment of occultation and the camera's spectral response (λ c ∼600 nm), the characteristic scale of the Fresnel diffraction is 1.2 km. To simplify the analysis, we ignore the Fresnel diffraction effect because it is so much smaller than our integration scale.

Description of Possible Detections
Six sites showed two or more consecutive flux drops of at least 2σ lower than the average light-curve value, all of those close to the nominal prediction time: (1-13) Bishop, (1-14) CPSLO, (2-06) Searchlight, (2-09) Mohave Valley, (V-07) Wildwood, and (2-11) Parker. These groups of flux drops can be seen in Figure 3 as two or more consecutive dropouts near relative time zero or to slightly negative or positive values inside the +1σ estimate. We note that several dropouts highlighted in the light curve of (2-14) Calipatria during the last 5 minutes are products of data taken through clouds, and those are not considered as possible detections. Upon closer examination, the flux drops from the six sites mentioned formed three separated groups when plotted in the sky plane. Figure 4 shows the area of interest around the possible detections in the sky plane, where each segment is a single integration and the segment color, from white to black, indicates the normalized flux from 0 to 1. The three groups of flux drops are labeled as 1, 2, and 3 in Figure 4 and in Figure 3. There is a group of flux drops seen in the light curves from Mohave Valley, Wildwood, and Parker about 30 s before the predicted time in Figure 3, labeled as 1, which translates in the group of flux drops seen in the sky plane around the coordinates (−600, −500)km in Figure 4. These stations are identified by red diamond symbols in Figure 1. A second group of flux drops is seen in the light curves from CPSLO, Searchlight, and Mohave Valley about 5 s after the predicted time in Figure 3, which translates in the group of flux drops seen in the sky plane around the coordinates (200, −300)km in Figure 4, and labeled as 2. These stations are identified by red diamond symbols in Figure 2. The data from Mohave Valley contain dropouts for both groups mentioned.
These two groups of flux drops, labeled 1 and 2, have a comparable maximum duration of ∼11 s (Searchlight and Wildwood sites), with a time elapsed between the two groups of flux drops of ∼36 s. These times translate to ∼250 km and ∼800 km, respectively, in the sky plane (see Figure 4). Additionally, the star does not disappear completely in the light curve from (V-07) Wildwood, (1-14) CPSLO, (2-06) Searchlight, and (2-09) Mohave Valley. Based on the S/N of the data (light curve and images), the residual flux at the bottom of the flux drops for those sites is inconsistent with zero, and the normalized flux drops to only about 50%. The similar duration and partial residual flux led us to consider models where these two groups of flux drops, labeled 1 and 2, correspond to an occultation of two similar-brightness stars with a projected separation of about 800 km at the distance of the occulting object (or about 38 mas). For comparison, we also consider alternative models where two objects occult a single star instead. For the same reason, this led us to not include in the comparison a model with a single object occulting a single star, as this cannot explain both groups of flux drops.
A third group of flux drops during at least two consecutive integrations was seen in the light curve from the (1-13) Bishop site. This translates in the sky plane to the flux drops seen around the coordinates (−600, −200)km in Figure 4, labeled 3. The Bishop site location is also indicated by a red diamond symbol in Figure 1. The shorter duration of these flux drops plus the additional constraint given by the nearby seemingly nondetections from (L-03) SwRI and (2-04) Indian Springs indicates a possible occultation of a secondary object. The light curve from Bishop is shown in more detail around the time of interest in Figure 5. The two upper panels show the raw light curves for two comparison stars; the third panel from the top shows the raw light curve for the target star. The lower panel shows the flux ratio between the target and the average of the comparison stars where this ratio has been normalized to unity for the unocculted flux. The predicted time is indicated by dotted-dashed green vertical lines, the 1σ in dashed green and . Light curves from all sites relative to the predicted time. The predicted time is indicated in dotted-dashed green vertical lines, the 1σ in dashed green and 2σ in solid green. Flux is normalized to unity when the star is not occulted. The right side shows the site code and name, cross-track predicted distance, and predicted time in UTC. An electronic copy of the data in this figure is provided. In red we highlight all data points higher or lower than 2σ, using the average S/N from each light curve. The labels 1, 2, and 3 identify the three groups of consecutive flux drops which are possible detections.
(The data used to create this figure are available.) 2σ in solid green. The midtime of the flux drops is about 30 s before the predicted occultation time, and it is separated in the sky plane from the other two groups of flux drops previously described by about 300 km and 750 km, respectively (see Figure 4).

The Occulted Star
Given the possibility that the occulted star could be a double (or an actual binary), we analyzed the publicly available astrometric data of the star. The best source of astrometry for this star is the astrometric stellar catalog of Gaia DR2. Unfortunately, multiple sources are not identified as such in Gaia DR2 and all sources are treated as single in the astrometric solution model (Brown et al. 2018). The so-called astrometric excess noise and the duplicated star flags, which could point to a multiple source, are both zero for this star. On the other hand, the empirical contrast sensitivity of Gaia DR2 derived by Brandeker & Cataldi (2019) indicates that a separation of 38 mas for equal brightness stars is about one order of magnitude below the capabilities of detection of multiple sources, so the apparent Gaia DR2 nondetection of this star as double is nonconstraining. There are other methods to detect binarity based on proper motion, but there are no data available to test them (Kervella et al. 2019). Solely from this information, the occurrence of a unresolved double star is plausible and cannot be discarded. Thus, we consider models with both cases, a single and a double star.

Occultation Modeling
We model the occulting object as a single object or as multiple objects, with all the objects in each model having either an elliptical or a circular projected shape. Given the presence of three groups of flux drops separated in the sky plane, for completeness, we consider models with up to three occulting objects. Each elliptical object has six free parameters: 1. the V-band absolute magnitude H V , 2. the geometric albedo in the V-band, p V , given by with au the astronomical units in kilometer, A the total projected area of one, two, or three objects depending of the model, and H e =−26.768 (Campins et al. 1985) the absolute magnitude of the Sun in the V band. We consider identical geometric albedo for all components in models with multiple objects. Together with H V , this uniquely defines the total projected area A of the object(s) and is the main result we seek from our analysis. 3. the center of the ellipse in the sky plane (x c , y c ). x c , y c are measured with respect to the nominal ephemerides of the object at the occultation time. 4. the minor semiaxis to major semiaxis ratio b/a of the ellipses, and 5. the ellipse position angle f of the major semiaxis measured from the north toward the east. For the simpler case of circular projected shapes, we have b/ a=1. For models with two and three objects, circular or elliptical, we consider identical geometric albedos for all object components, and we introduce an additional parameter that defines the fraction of the total projected area A occupied by the primary A p , secondary A s , and tertiary A t components with the constraints for two-and three-object models, respectively. As mentioned before, for the models with multiple objects, we consider only the case where all the components have either circular or elliptical projected shapes, and we do not include combinations of circular and elliptical components. This would increase the number of models to compare, and ultimately, we are mostly interested in analyzing whether the data favor a model with single or multiple objects, and whether the data favors a more complex model with elliptical projected shapes versus simpler models with circular projected shapes.
For the star, we consider the occultation of one or two stars. With this option, our model will have the freedom to adapt to the nonzero flux during an occultation. The stars are modeled as point sources where the normalized flux is unity when not  ( Model with no occultation .. 0 A model where the star is not occulted and it is assumed that the variations in the light curves are due entirely to random Gaussian noise. The model light curves are flat with normalized flux equal to one at all times for all the sites.
Models where only two groups of flux drops, labeled 1 and 2, are produced by an occultation A single object with a circular projected shape occulting two stars that models two groups of flux drops. The flux drops seen in the Bishop site are considered to be just due to noise and not a real occultation.
Two objects with circular projected shapes occulting a single star that models two groups of flux drops. As with  1 , the flux drops seen in the Bishop site are considered to be just due to noise and not a real occultation.
Similar to  1 but here the object has an elliptical projected shape.
Similar to  2 but both objects have elliptical projected shapes.
Models where the three groups of flux drops are produced by an occultation Two objects with circular projected shapes occulting two stars. The detection from Bishop is now modeled as a secondary object occulting the primary star.  6 2 (both circular) Two objects with circular profiles occulting two stars. Similar to  5 but provides an alternative geometry where the detection from Bishop corresponds to a secondary object occulting the secondary star.
Three objects with circular projected shapes occulting a single star that attempts to match the three groups of flux drops.  8 2 (both elliptical) Similar to  5 but both objects have elliptical projected shapes.  9 2 (both elliptical) Similar to  7 but the three objects have elliptical projected shapes.
Note.The first column is the model identification. The second column is the number of occulting objects modeled (0, 1, or 2). The third column is the list of free parameters for the occulting object. The subscripts "p," "s," and "t" refer to the primary, secondary, and the tertiary components, respectively. The fourth column is the number of occulted stars modeled (1 or 2). The fifth column is the list of free parameters for the occulted star(s). The subscripts "p" and "s" refer to the primary and secondary stellar components, respectively. The sixth column is the total number of free parameters, which indicates the model complexity. The seventh column is a description of the model. Models are grouped such that the model where there is no occultation and all the flux drops are produced by random noise is in the upper section. In the central section are the models where the two groups of flux drops seen in Figure 3, labeled 1 and 2 in Figures 3 and 4, are the product of a real occultation. In the lower section are the models where all three groups of flux drops are product of a real occultation, including the flux drops seen in the Bishop site. In each section, the models are ordered by the number of free parameters.
occulted. For models with two stars, we introduce three additional free parameters: 1. the relative position of the secondary star with respect to the primary projected at the occulting object distance (x star , y star ) s , and 2. the flux of the primary star f p .
We note that, given the faintness of the occulting object, the flux of the secondary star ( f s ) is constrained to be f s =1−f p .

Alternative Models
The possibility of one, two, or three objects occulting one or two stars leads to 10 different models (labeled  1 to  10 ) to compare, each with a different complexity driven by the number of free parameters. For completeness, we add to the model comparison a model  0 where there is no occultation. In model  , 0 all of the flux drops are the product of random noise fluctuations in the light curves. The models vary in complexity from 7 up to 16 free parameters while the model  0 has no free parameters. These models are explained in detail in Table 4.
The first column in Table 4 is the model identification. The second column is the number of occulting objects in the model with the type of projected shape indicated in parentheses. The third column lists the free parameters for the occulting object (s). The subscripts "p," "s," and "t" refer to the primary, secondary, and tertiary components, respectively. The fourth column contains the number of occulted stars in the model. The fifth column lists the free parameters for the occulted star(s). The subscripts "p" and "s" refer to the primary and secondary stellar components, respectively. The sixth column provides the total number of free parameters in the complete model. The seventh column is a description of the model.

Statistical Analysis
The traditional approach to analyze stellar occultation by opaque objects (without atmosphere) has been (1) to model each occultation light curve with a square model including, when appropriate, the effects of the star angular size and diffraction, (2) determine the times of ingress and egress of the star behind the object as seen from each observing station by fitting a square model, and (3) to compare the times of ingress and egress measured with the expected ingress and egress times for a model of the object. This relatively simple approach is sufficient for occultations where the occultation is unambiguous because of light curves with a high S/N, long duration of the occultation with respect to the exposure time, or both. The approach becomes less practical in cases where the presence of an occultation is the ambiguous product of a short event duration with respect to the exposure time, low S/N, or both. In this work, we adopt an approach where the measured flux in each data point of the light curves is compared with the flux of a model light curve, avoiding the need to fit against formal times of ingress and egress. Consequently, there is no need to define a priori which data point is or is not part of an occultation, which is analyzed instead as part of a Bayesian inference framework. For instance, these ingress and egress times may not exist at all if the putative detections are produced by random noise instead of an occultation by the object. The adopted approach naturally covers cases of seemingly negative detections where the object is grazing and the magnitude of the flux drop in the light curve is below the noise level.

Bayesian Inference
Given the possibility of different models, we require a quantitative way to compare them and to determine which one is statistically favored given the occultation data, while taking into account the model complexity. In this work, we adopt a traditional two-step Bayesian inference approach (Gregory 2005).
Step one of the inference consists in comparing and quantifying the support for each model over another given the occultation data and prior information. This step uses the marginal likelihood   D k ( | ) of each model, also known as evidence, model evidence, and global likelihood, to quantitatively rank the alternative models and to determine the favored model, the one with the highest marginal likelihood. The marginal likelihood   D k ( | ) is obtained by integrating the likelihood function weighted by the prior of the parameters over the whole parameter space of the model (Gregory 2005): where the data D are the normalized fluxes f i as a function of the midtime t i of each integration of the light curves, θ are the parameters for the model  k , q   D , k ( | ) is the likelihood function, and P(θ) is a prior distribution for the parameters θ. We note that the marginal likelihood can only be used to compare whole models when the same data are considered in each model (Edwards 1992). As such, the marginal likelihood   D k ( | ) is a number that can be used to rank the alternative models, with the model with the highest   D k ( | ) being favored over the alternative ones. We recall that, from the Bayes rule, the posterior probability  P D k ( | ) of a model  k given the data D is given by (Gregory 2005 ) is the prior probability for the model, and P D ( ) is a normalization constant. From this, any pair of models  k and  m can then be statistically compared with the use of the posterior odds, given by Given that there is no observational data on the rate of close and small binaries or on the incidence of double stars with a small angular separation, we assume that all alternative models Note.The first column is the range of values for the Bayes factor, and the second column is the empirical interpretation of the value from the Jeffreys scale. Reproduced from Jeffreys (1998).
are equally probable a priori, in which case we have and the models can be statistically compared with the Bayes factor  km , given by the ratio of the marginal likelihoods We adopt the empirical Jeffreys scale (Jeffreys 1998) to assist in the interpretation of the Bayes factor between pairs of models, reproduced for convenience in Table 5. The likelihood function for normally distributed uncorrelated flux uncertainties is given by (Gregory 2005) is the modeled light curve given the parameter vector θ of the model  k , σ i is the uncertainty in the flux of the ith data point, and N is the total number of data points from all the light curves. The uncertainties in the flux σ i are derived from the photometry and modeled as a zero-mean normal distribution. The total number of data points N used for the analysis is given by where n j is the number of data points used from the jth light curve corresponding to the jth site, and N LC is the total number of light curves used.
We numerically estimate the marginal likelihood   D k ( | ) for each model using the python package dynesty (Speagle 2020), an implementation of the nested sampling algorithm (Skilling 2004(Skilling , 2006. For practical reasons, dynesty returns a numerical estimation of the natural logarithm of   D k ( | ), which is the value that we report below in our results (see Table 6). The sampling with dynesty is done with 512 live points, which gives robust estimations of   D k ( | ), i.e., a larger number of live points does not change the estimated value of   D k ( | ) nor their formal uncertainties. For the bound method in dynesty, we adopt multiple (possibly overlapping) bounding ellipsoids. The dynesty package provides other bounding methods, although we found no systematic difference in the results or performance when a single bounding ellipsoid is used. For each sampling, we initialize the live-point positions distributed uniformly in the parameter space defined by the priors. Details about the nested sampling algorithm, the dynesty implementation of this algorithm, and the meaning of the live points can be found in Skilling (2004Skilling ( , 2006 and , respectively. Step two of the Bayesian inference consists in estimating the posterior probability density function (pdf) q  P D, k ( | ) of the parameters θ of the favored model  k . From the posterior pdf, we obtain the nominal parameter values and the formal uncertainties for the favored model. We recall that, from the Bayes rule, we have (Gregory 2005)  Note.Results of the model comparison, in the same order shown in Table 4. The statistically favored model ( 5 ) is highlighted in bold letters. The tabulated Bayes factor  k5 for other models is referred to this. The first column is the model identification from Table 4. The second column is the number of occulting objects modeled (0, 1, or 2). The third column is the number of occulted stars modeled (1 or 2). The fourth column is the total number of free parameters. The fifth column is the natural logarithm of the marginal likelihood   D ln k ( ( | )) estimated with the nested sampling algorithm, except for the model without occultation  0 . The sixth column is the Bayes factor between each model  k and the favored model  5 . The last column is the empirical interpretation adopting the Jeffrey's scale (Jeffreys 1998). where q   D , k ( | ) is the likelihood function, q P ( ) is a prior distribution for the parameters θ, and   D k ( | ) is the marginal likelihood.
For the numerical estimation of q  P D, k ( | ), we use a Markov Chain Monte Carlo algorithm to obtain representative samples of the posterior pdf. We adopt the use of the Python package emcee 3.0 (Foreman-Mackey et al. 2013), which implements the affine-invariant ensemble sampler (Goodman & Weare 2010). For each model, we use the sampler with n w =512 random walkers with initial positions distributed uniformly in the parameter space defined by the priors. The implementation in emcee of the integration autocorrelation time t iac is used to evaluate the chain convergence and to determine the length of the so-called burn-in stage, measured in number of steps. The burn-in stage consists of 500 steps, which comply with the previous criteria. After this, the chain is run for an additional n step =200 steps to obtain n step ×n walker =102,400 samples. The samples from the burn-in stage are discarded, and the joint and marginal probability distributions for the parameters are estimated from the histograms of the remaining samples. The nominal parameter values of the model are obtained from the maximum of the marginal posterior pdf of the parameters of interest. The formal uncertainties are taken from the 68% credible intervals in the marginal posterior pdf.

Parameter Priors
We specify the prior information for the model parameters P (θ) using probability distributions. For the absolute V-band magnitude H V , we adopt a normal distribution with mean of 7.2 and standard deviation of 0.3 estimated from the photometry retrieved from MPC. We do not attempt to account for the inevitable systematic errors in such photometry but leave that for a later improvement in a follow-on project if deemed important enough. For the geometric albedo in the V-band p V , we adopt a normal distribution with mean value μ=7.2%, standard deviation σ=7.6%, and truncated between 2% and 100%, fitted to the empirical albedo distribution of TNOs from Kovalenko et al. (2017). Our use of H V and p V means that the total projected area A of the occulting object(s) is not a direct free parameter in the model and is instead computed from the albedo and absolute magnitude, using the relation in Equation (1).
For the center of each object (x c , y c ) in the sky plane, the prior distributions are chosen to be uninformative while isolating each of the flux drop groups (see Figure 4) and to maintain the identity of each object component during the sampling of the parameter space. We choose uniform distributions constraining the center positions to a box 200 km around the approximate center of each possible detection, as seen in the sky plane (Figure 4). The size of this box is the same for all the models to enable a proper model comparison. For the fractional areas A p , A s , and A t , we adopt uninformative priors with a uniform distribution between 0.0 and 1.0. Recall that with multiple objects, the total fractional area is required to sum to unity. We note that for models with a single object, this constraint reduces to the trivial constraint of A p =1. The priors for the orientation of the ellipses, f, are uninformative uniform distributions between 0°and 180°.
For models with two stars, we define a prior for f p , which constrains the relative fluxes of the stars. Given that we have no available constraints on the relative fluxes of the two stars, the prior for f p is chosen to be an uninformative uniform distribution between 0 and 1. We note that for models with a single star, this reduces to f p =1, meaning that the normalized flux goes to 0 when the star is occulted, recalling that this is due to the fact that the apparent magnitude of the object is much fainter than the occulted star and the magnitude limit of the images (V∼17), making the object flux well below the background noise level during an occultation. For the relative position of the secondary star with respect to the primary, we chose a uniform distribution constraining the center position of the secondary to a box 200 km around the approximate center of the flux drops, as seen in the sky plane (Figure 4). Table 6 summarizes the result of the statistical comparison of the models. For clarity, the first four columns repeat some information of the alternative models from Table 4: the model identification, the number of occulting objects in the model, the number of occulted stars in the model, and the number of free parameters, respectively. The fifth column informs the natural logarithm of the marginal likelihood   D ln k ( ( | )) estimated with the nested sampling algorithm, except for the model without occultation  0 that is calculated directly from Equation (9). From this analysis, the conclusion is that the model  5 is statistically favored given the occultation data. Thus, we infer that our observation is of a binary TNO occulting a binary (or double) star. For clarity, we choose to give in the sixth column only the Bayes factor  k5 between each model  k and the favored model  5 . Any other pair of models k and m can be compared with Bayes factors using either the marginal likelihood values, or, for instance

Model Comparison
) and Bayes factors in Table 6, we can draw additional inferences: 1.  0 with no occultation can be discarded, and there is decisive evidence for an occultation present in the data. This is not surprising but nonetheless reassuring. 2. There is decisive evidence against all the alternative models with only one occulted star ( 2 ,  4 ,  7 ,  10 ). This is not a surprising result due to the depth of the flux drops seen in the light curves. 3. The choice between  5 and  6 comes down to which stellar component is occulted by the secondary body as shown in Figure 6. The favored model requires a nondetection of the secondary body occulting the secondary star at SwRI and Bishop (see the upper panel of Figure 6). The nonfavored model requires a nondetection of the secondary body occulting the primary star at Bishop and Indian Springs (lower panel of Figure 6). As seen in the figure, the latter case requires a much tighter fit of the object between the tracks and is correspondingly less likely compared to the favored case. Our inference in this case is still substantial but not quite as strong. When using these results for future occultation predictions, it may be wise to consider both alternatives geometries of models  5 and  6 weighted by their respective marginal likelihoods. 4. There is only weak evidence against the model  8 , which is similar to the favored one, but with the two objects having elliptical projected shapes. Our interpretation is that the additional model complexity due to the extra free parameters is not supported by the data. The penalty incurred from the additional free parameters is larger than any improvement in the capacity of the model  8 to predict the occultation data, given the data S/N. Still,  8 is the second-most favorable model. Higher S/ N data from upcoming occultations or rotational light curves could be used in the future to bring new and stronger evidence about the projected shape and tridimensional shape of the objects.

Physical Properties of the Occulting Object
We proceed in turn to estimate the parameter values and formal uncertainties for the statistically favored model ( 5 ) identified in Section 8. Figure 7 shows the posterior pdf for the geometric albedo and the diameter for each component indicating that the binary object is dark with a geometric albedo p V =5.1%±1.7%, a primary of diameter D=181±16 km, and a secondary of diameter D=138±32 km. Figure 8 shows the posterior pdf for the location of the secondary component with respect to the primary. The separation is For illustration, Figure 9 compares the normalized light curves with modeled light curves calculated using the nominal parameters values for the model  5 ( Table 7). The modeled light curve are shown by the empty green square symbols while the occultation data with their formal uncertainties are shown in black circular symbols. For each site, the normalized flux is plotted as a function of time referred to the predicted time. The nominal light-curve model indicates that the primary component occulting the primary star is detected in the (2-09) Mohave Valley, (V-07) Wildwood, and (2-11) Parker sites while the primary component occulting the secondary star is detected in the (1-14) CPSLO, (2-06) Searchlight, and (2-09) Mohave Valley sites. The secondary component occulting the primary star is detected only in the (1-13) Bishop site. We note that the illustrated light curve is for the nominal parameter values of model  5 , but the model  5 as a whole involves a probability distribution for the parameters instead, which is an advantage of performing a Bayesian analysis of the occultation data. For the same reason, we cannot indicate which data point is part of the occultation, because each data point has a    probability of being part of the occultation, the product of the model parameters being a pdf and not a single value. The top panel of Figure 6 shows the nominal geometry, using the nominal parameter values, for the statistically favored model  5 projected in the sky plane, with the x-axis in the east direction and the y-axis in the north direction. The primary and secondary objects occulting the primary star are drawn in solid green circles in the lower right and upper right, respectively. The dashed blue circles show the same object but occulting the secondary star. Object radii and separation are with the nominal parameter values from Table 7. Figure 1 shows the ground tracks of the object occulting the primary star using the nominal values from Table 7. Figure 2 shows the ground track of the objects occulting the secondary star. The stations with gray squares are those with no data or unconstraining data (late start or wrong field), in blue circles those stations with no detection (with the nominal parameter values of the favored model), and in red diamonds those stations with a detection. Figure 10 shows 2014 WC 510 compared with other TNOs (excluding Centaurs) with published and well-determined equivalent diameter d eq and geometric albedo p V (relative uncertainties smaller than 10% in effective diameter). The primary and secondary components are indicated by black squared symbols labeled P and S, respectively. The larger black circular symbol indicates the equivalent diameter considering both components. The remaining solid symbols are measurements from published stellar occultations:  Arrokoth (Buie et al. 2020a), Leleākūhonua (Buie et al. 2020b), while empty symbols are those from other methods (mainly thermal) selected from the compilation by Johnston (2018). Besides Arrokoth and Leleākūhonua, 2014 WC 510 is the smallest object measured from a stellar occultation, similar in size to Leleākūhonua and similar in size and albedo to the resonant object 2003UT 292 measured with radiometry. Figure 11 shows the posterior pdf for the geometry of the occulted stars. From top to bottom, the diagonal panels are the marginal posterior pdf for the angular separation, position angle, and flux ratio of the secondary star with respect to the primary star. We note that the a priori choice to dub each star primary and secondary was arbitrary and from the flux ratio in the posterior pdf of Figure 11, the secondary could be actually brighter than the primary. The angular separation of the stars is well constrained to 38±1 mas driven by several chords detecting the primary object. An independent confirmation of the duplicity of the stars' separation and measurement of the flux ratio is within the theoretical capabilities, although in the limit of speckle imaging from an 8 m class telescope, but at the time of this work, these data are not available.

Astrometric Constraints
A valuable by-product of a stellar occultation is that the occulting object has a known position with respect to the occulted star at the time of the occultation with an accuracy no worse than the object size and often much better. In turn, this astrometric position can be used to improve the object's orbit fit and, consequently, reduce the uncertainty of future stellar occultations. Here we analyze this astrometric constraint discussing the systematic uncertainties introduced by the duplicity of the occulted star. It is worth mentioning that the geometry of the double star components derived in this work is at the time of the occultation and differs from the Gaia DR2 stellar catalog epoch by 3.4 yr. At this moment it is not possible to know if the stars' geometry was different and by how much when it was measured by Gaia. Similarly, any measurement from speckle imaging from the ground will have a different epoch from the occultation, and the double star could have a different geometry by then, either facilitating or complicating an independent confirmation of its double nature.
The favored model  5 considers two stars, and we adopt that the astrometric position of the star from Gaia DR2 ( Table 1) is the photocenter of the two stars. Figure 12 shows the posterior pdf for the position of the primary star and secondary stars with respect to the star photocenter. The photocenter is at 17.5±2.2 mas from the primary star. From this, we derive the astrometric position of each star at the occultation time ( Table 8). The uncertainties in Table 8 do not include the uncertainties of the star from Gaia DR2 given in Table 1.
A source of systematic error in the position of the double star components is the unknown nature of the occulted stars. The stars could be an actual bounded binary system or a visual double. If the star is an actual bound binary, then the Gaia DR2 proper motion is a combined effect of the proper motion of the star system barycenter and the relative movement of the photocenter around the barycenter. On the other hand, if the stars are not bound and instead aligned by chance, then the  Gaia DR2 proper motion is a combined effect of the proper motion and parallax of each star. Treating the stars as bound, the angular distance of 38 mas derived from the occultation analysis translates in a physical projected distance of 68 au, which is well in the range of known binary stellar systems. Nonetheless, it is not possible at this moment to discard the possibility of a chance alignment of two unbounded stars at different distances. Discriminating between these two possibilities goes beyond the scope of this work but could be resolved with direct imaging measurements of the stars, if that is possible.
With the adopted astrometric position of the stars, we calculate the position of the occulting body's primary and secondary components at the reference time t 0 =2018 December 1 12:21:23 UTC and given in Table 8. The uncertainties in the astrometric position of the object components include the uncertainties in the double star parameters and from the stellar catalog at the occultation epoch from Table 1.

Discussion
The statistical analysis of the data from the first stellar occultation by the resonant object 2014 WC 510 favors a model of a binary object with a primary of diameter d p =181 km and a secondary with diameter d s =138 km. The projected shape of the object components is circular with weak statistical evidence against elliptical projected shapes, shapes that could be confirmed with higher S/N data in future occultations, or rotational light-curve data. This object is among the darkest objects measured with stellar occultations with a geometric albedo in the V band of 5%. This albedo is clearly lower than for the cold-classical Arrokoth and happens to be comparable to values typical for the Jupiter Trojans. This albedo measurement could be subject to a systematic error depending on its light-curve properties. We only know the mean brightness and the projected circular shape derived from the occultation. If the tridimensional shape of the object components are spherical, then our measurement should be accurate. If there is a rotational or a secular light curve from nonsphericity or from phase effects, then there will be a need for future recalibration of the albedo once better photometry is available. We note that the occultation data favor a simpler model with circular projected shapes over a more complex model with elliptical projected shapes. These circular projected shapes could be a product of spherical object components or the product of elongated object components, which can give circular or close-to-circular projected shapes depending of the  Table 8. The photocenter is the astrometric position derived from Gaia DR2 (see Table 1). The contours in the joint pdf's are the 39.3% credible intervals. (The data used to create this figure are available.) Note.Astrometric constraints derived from the occultation using the preferred model M 5 . The star component positions and uncertainties are derived from Figure 12. All astrometric positions are given in the J2000 reference frame. The astrometric position of the object components are for the reference time t 0 =2018 December 1 12:21:23 UTC.
orientation of the rotation pole and rotation phase, possibilities that the present data cannot discriminate. From the statistical analysis of the occultation data, the partial disappearance of the star favors a model where the occulted stars is a previously unresolved double star (a visual double or an actual binary) with similar-brightness components and separated by 38 mas. With the available data, it is not possible to determine how different was the geometry of the stars when measured by Gaia DR2, which introduces a systematic uncertainty in the position of the star components at the occultation time, which in turn propagates to the derived position of the occulting object. More accurate parameters for the star could be available in future releases from the Gaia mission or from ground-based follow-up observations, such as speckle imaging to resolve both components, although 38 mas is in the limit of the theoretical capabilities of speckle imaging for a 8 m class telescope (Matson et al. 2019). Speckle imaging could provide, in principle, an independent measurement of the stars' relative brightness and separation to improve the modeling of the double object components from the occultation. However, it must be taken into account that attempts to confirm the double nature of the occulted star will also suffer from the unknown nature of the stars (double or actual binary), and the projected separation and position angle of the star components could also be different with respect to the one derived in the present occultation, when and if finally measured. To improve the characterization of 2014 WC 510 , we consider it more valuable to pursue future occultations together with studies of its rotational light curve. Observers must take into account the possibility of double stars present in the Gaia DR2 catalog not being labeled as such. The empirical contrast sensitivity of Gaia DR2 (Brandeker & Cataldi 2019) gives good constraints in the possibility of these unresolved double (or multiple) stars to take into account in the planning of occultation campaigns, which could be considered as systematic biases in the predictions.
This occultation measurement puts 2014 WC 510 at the small end of the known multiple TNOs with accurate size and albedo measurements, with the exception of the small classical TNO Arrokoth (Stern et al. 2019;Buie et al. 2020b). Its double nature is not surprising and supports the expectation that most TNOs were formed as binaries, and that closer binaries have more chances of survival (Fraser et al. 2017;Nesvorný & Vokrouhlický 2019). The derived diameter of the object components put 2014 WC 510 in the category of similar-size binaries, which could have a primordial origin, unlike the large TNOs with small satellites thought to have formed after catastrophic impacts (Noll et al. 2020). Stellar occultations by TNOs is a promising technique to advance the characterization of binary TNOs, particularly for the tight and small systems which are inaccessible by direct imaging.
2014 WC 510 is moving in front of the galactic plane during the next few years, providing occultation opportunities to refine the size, shape, and orbit. From the RECON event prediction list, there are 30 occultation opportunities coming up in the next seven years visible somewhere on Earth, with seven of these visible from RECON itself. In future efforts, we recommend station separations smaller than 50 km to sample the secondary object. The overall coverage needed to ensure capturing both objects is not known at this time. The positions measured here can only provide guidelines, but it is clear the coverage required could still be as large as 1000 km due to the eccentricity of the mutual orbit even if the center of mass can be predicted with perfect accuracy.
This measurement of 2014 WC 510 helps extend the coverage in the size-albedo parameter space in which good measurements are absent below ∼200 km, with the unique exception of Arrokoth. The characterization of small objects in this range (<200 km) is essential to understand the composition and evolution of the small population of TNOs. 2014 WC 510 is consistent with the general tendency of smaller objects to have darker surfaces. A larger sample of accurate albedo and size measurements is necessary to determine if this tendency is real or due to an observational bias effect against higher albedo objects.