The Complex Rotational Light Curve of (385446) Manwë–Thorondor, a Multicomponent Eclipsing System in the Kuiper Belt

The epoch 1 observations reported here were obtained in the Sloan Digital Sky Survey (SDSS) r and g bands on two consecutive nights (2016 August 25 and 26 UT) using the SOAR Optical Imager on the 4 m SOAR telescope at Cerro Pachón, Chile. The imager consists of a minimosaic with two thinned, back-illuminated E2V 2k × 4k charge-coupled devices (CCDs), together covering a 526 square field. We used the 2 × 2 binned readout mode, yielding a pixel scale of 0154. On both nights we took advantage of the remote-observing capabilities of the SOAR telescope, controlling the camera and telescope from our home institutions in the USA. Operators and staff at the observatory provided instrument support and monitored the telescope, camera, and weather status throughout the night.

On both nights we observed Manwë–Thorondor at 8–10 minute intervals over a period of 8 hr (02:00–10:00 hr UT), with a nominal observing cadence of four r-band exposures followed by one in g band. Most exposures were fixed at 450 and 600 s for the r and g bands, respectively. Except for small dithers between exposures, we targeted the same field on both nights. This allowed us to use an identical set of field stars to calibrate all the photometry. The two-night motion of Manwë–Thorondor was 15, well within the field covered by the camera. We also oriented the camera so that the target did not move across the gap between the two CCDs in the array. Both nights were clear, with seeing variable from 05 to 20 (median 12).

Our second night was nearly optimally timed to cover the mutual event with an expected midtime, t_mid = 05:25 hr UT August 26, and duration ∼10 hr. Here we define mutual event duration to be the interval of time when the expected brightness diminution would exceed 0.1 mag (see Grundy et al. 2014). Given the 12 hr uncertainty in the timing and our ∼8 hr observation windows, our likelihood of detecting at least half the event was ∼29% (5.3% and 24.1% on the nights of August 25 and 26, respectively). For detecting at least a 10% diminution, our expected likelihood was ∼64% (13% and 51% on the respective nights).

The epoch 2 observations were obtained in SDSS r and g bands on 2017 July 22 and 25 UT using the GMOS-S imager on the 8.1 m Gemini South telescope at Cerro Pachón, Chile (the second night was delayed by 2 nights due to bad weather). The time was awarded under Gemini's "Fast Turnaround Time" program. The imager is an array of three Hamumatsu 2048 × 4176 pixel CCDs with a field of 6144 × 4224 pixels covering 8.35 × 563. With 2 × 2 binning, the pixel scale is 016. The queued observations were obtained by onsite operators on both nights and sequenced continuously while conditions were acceptable (most of night 1, part of night 2). We monitored the operations remotely with a Skype connection.

Sequences consisted of five r-band exposures (240 s each) followed by two g-band exposures (150 s each), repeated throughout the night. All exposures were dithered, with the same pointing used for all exposures on both nights. The target was centered in the field on the first night but appeared near a gap between two detectors on the second night. Because of the dithering on the 2nd night, the target sometimes fell into the gap between two of the detectors. Both nights were marginally clear with seeing variable from 10 to 20 and with median values 12 (night 1) and 16 (night 2).

The epoch 3 observations were obtained on 2017 November 9 UT in r band using the Large Monolithic Imager (LMI) on Lowell Observatory's Discovery Channel Telescope (DCT) at Happy Jack AZ. The LMI is a single 6144 × 6160 pixel E2V CCD with 024 pixel scale (3 × 3 binned) covering a field 12.3 × 123. Observations had been planned for the latter halves of both 2017 November 8 and 9 UT. Because of bad weather, only sparse observations were obtained, and only on the second night. All exposures were 600 s. Seeing was 1'' to 15.

The epoch 4 observations were obtained 2018 October 6 UT in a wide VR band, again using the LMI on the DCT. Seeing was stable, between 12 and 15. The target was too faint to observe using narrower-band filters.

To obtain the SOAR and Gemini light curves, we first used standard techniques to remove bias and flat-field artifacts from the raw images, relying on bias exposures and twilight flats taken on both nights. For both the SOAR and Gemini data sets, each consisting of two nights of observations, we then used two different procedures to calibrate the data. With procedure 1 we used IRAF programs to obtain a relative calibration of the target magnitudes, with careful attention to account not only for the relative extinction and seeing from one exposure to the next but also for the influence of faint field stars near the path of the target. With procedure 2 we obtained absolutely calibrated magnitudes, including color corrections, making use of a photometry pipeline that measures and calibrates source magnitudes for fields that are within the footprint of the SDSS. This second procedure is described in Rabinowitz et al. (2014).

We describe our implementation of procedure 1 in detail below.

Step 1: We determined the magnitude zero-points for all the images on both nights relative to a chosen best image (least extinction, best seeing) by measuring the instrumental magnitudes for a fixed set of field stars common to all the images. We used a large aperture with diameter nearly double the worst seeing (25 for SOAR, 32 for Gemini) for the field star photometry, and also a smaller aperture comparable to the median seeing (12 for SOAR, 13 for Gemini) to determine an aperture correction for the target photometry (see step 4).

Step 2: For each night we constructed a median reference image by selecting the best images (least extinction, best seeing, deepest exposures), sky-subtracting them, shifting their pixels to achieve a common alignment, scaling them to compensate for their relative zero-points, and then computing the median average at each pixel location.

Step 3: For each image acquired on a given night, we removed the field stars by subtracting the alternate night's median reference image after shifting each image to align with the reference, and scaling the reference to have the same zero-point as the image. An alternate night's reference image was chosen to make sure there was no residual signal in the reference image at the location of the moving target.

Step 4: Using the same small aperture used in Step 1 to measure the field stars, we measured the magnitudes for the target in the reference-subtracted images. We then used relative zero-points and aperture corrections determined from the field stars (step 1) to normalize the photometry.

For the r band, we used procedure 1 to obtain a precise measure of the light curve shape and then relied on procedure 2 only to have an independent check on these measurements and to establish the magnitude offset and color terms required to transform the photometry into the Sloan system. For the g band, because we had only one night of good data (Gemini) or because we did not have enough observations to construct a good median sky image for both nights (SOAR), we relied exclusively on the second procedure for all photometry.

For both procedures 1 and 2 we used a small aperture comparable to the median seeing to measure the target flux. Note that the motion of the target in the SOAR r-band images (nearly all with 450 s exposure times) was 034, leading to a negligible influence on the photometry (<1% loss in signal) given the much larger photometric aperture (15). The g-band exposures were longer (600 s), leading to a possible small trailing loss of ∼1.5%. There was insignificant trailing loss in the short Gemini exposures.

For both the r and g photometry, we determined color corrections iteratively in two passes under the assumption that the corrections were the same for all images in a given band. In the first pass we made an initial estimate of g–r without color corrections, knowing this would require a later correction. In the second pass we used our estimate for g–r to determine the color corrections individually for each measurement. Then for each band we evaluated the mean color correction and used this to correct all the measurement from both nights. Our final estimate for g–r is derived from these color-corrected values for g and r.

The DCT observations were reduced using the standard procedures described by Thirouin & Sheppard (2019).

The r-band data collected on 2018 November 8 were taken through clouds in poor conditions. We include the few observations that could be calibrated to standards but give them no weight in the foregoing analysis.

The VR-band data were calibrated to r-band standards. In the analysis below, we attempt to simultaneously fit the resulting r-band light curve along with our r-band photometry from previous epochs to a common model, ignoring any possible changes in color over time that might affect the relative calibration of the r and VR bands to the same magnitude system.

Table 3 lists all our magnitude measurements, M, along with Julian Date (JD) at the start of each exposure, the measurement error, the exposure time of the associated image, and the telescope. Light travel-time corrections relative to a common epoch (JD 2454300.0) were computed for each observation date and used in the analysis described below. For the r band, the listed measurement errors do not include the systematic uncertainty from the absolute zero-point calibration and color correction. They do include the uncertainties in the instrumental magnitudes and their relative normalization determined from the statistical noise in the measured sky and target flux. The zero-point uncertainties for both the g- and r-band magnitudes are ∼2%, based on unsaturated stars in the field with SDSS magnitude uncertainties <10%. Measurements contaminated by cosmic ray hits or taken in poor observing conditions (wind bounce, high airmass, poor focus, or insufficient exposure time) are excluded from the tables and the analysis described below.

Table 4 lists resolved photometry of Manwë and Thorondor (magnitudes M_M and M_T, respectively) previously obtained with the HST Wide-Field Planetary Camera2 (WFPC2) in the F606W band at seven different epochs from 2007 to 2013 (earlier observations on 2006 July 25 are excluded because no F606W-band data were obtained). Given the allocated orbits, only a few observations were obtained at each epoch, with each epoch spanning <0.1 day. Grundy et al. (2014) report magnitudes for Manwë and Thorondor derived from these data after transforming F606W to V band and averaging the results from each epoch. Here we present the measurements without transformation or averaging. For each epoch, we also list θ and the apparent separation, Δ, between Manwë and Thorondor. The listed magnitude errors are the variance with respect to the epoch mean. For the 2008 September 7 epoch, there is only a single measurement. We estimate an error of 0.1 mag consistent with the error for other epochs. An asterisk (*) indicates rejected measurements in the foregoing analysis (outliers, measurements with unusually large variance, or measurements made when Manwë and Thorondor were too close together to deconvolve reliably).

For each HST observation, we also list a pseudo magnitude for Manwë+Thorondor, M*, which is the F606W magnitude that would be observed for Manwë+Thorondor (unresolved) if the magnitude of Thorondor were fixed at its weighted mean value, $\langle {M}_{{\rm{T}}}\rangle =24.4$ mag, determined from all the unrejected HST observations. The reason for constructing this derived signal will be discussed further below. Here we note that if Manwë's light curve is regular and periodic on a timescale of ∼1 day, while Thorondor's brightness varies on a timescale comparable to the time between observing epochs (months to years), then a single rotational model will fit both the long-term variability of M* observed with HST in the years 2007 to 2013 and the more recently observed short-term variability we observe for the unresolved system with the SOAR, Gemini, and DCT telescopes.

Here we determine the rotation period P by simultaneously fitting a periodic function, M(t), to all the photometry listed in Table 3 with the exception of the first DCT epoch (2017 November 8) which were collected in bad observing conditions, as noted above. We also fit the unrejected combined HST photometry of Manwë+Thorondor listed in Table 4, M*. As discussed above, we construct M* assuming a brightness for Thorondor fixed at its mean HST brightness. Our assumption, here, is that the light curve of Manwë is regular and periodic over the entire time period spanned by the observations, while the variability of Thorondor is on a timescale much longer than any given observation epoch. The constructed HST photometry would therefore have the same periodicity and shape as the average light curve observed over a single epoch from the ground.

Following Rabinowitz et al. (2014), we express M(t) as

$$\begin{eqnarray}&&M(t)={H}_{{\rm{j}}}+\beta \theta (t)+F(t)\end{eqnarray} \tag{ 1 }$$

where H_j is the mean reduced magnitude in band j at phase angle θ = 0 (i.e., the absolute magnitude in each band), t is the light-travel corrected observation time, and β parameterizes a linear dependence on θ. Function F(t) is a single-harmonic sinusoid with harmonic amplitude A and phase ϕ,

$$\begin{eqnarray}&&F(t)=(A/2)\sin (\pi [2\phi +4\omega (P,t)])\end{eqnarray} \tag{ 2 }$$

where ω(P, t) is the rotational phase, ω(P, t) = (t–t_o)/P, referred to an arbitrary epoch, t_o. There are two peaks per rotation, corresponding to the light curve of a rotating ellipsoid.

The first row of Table 5 lists the minimum χ², the number of degrees of freedom, N (239 observations less 7 fitted parameters), and the best-fit values and uncertainties for P, β, A, ϕ, the absolute magnitudes H_g, H_r, and H_F606W, and g–r. Here, we have assumed t_o = JD 2544300.0, and we use the χ² minimization procedure described in detail by Rabinowitz et al. (2014). We also used Monte Carlo simulations to determine the 68% confidence intervals of the associated parameters. These simulations preserve the observation times and measurement errors but use the best-fit light curve to model the magnitude values. These confidence intervals are valid if the best-fit light curve is a good representation of the observed signal (i.e., χ² ∼ 1). The second row of parameter values in Table 5 refers to an improved fit after making additional adjustments to the data. We discuss the differences in the resulting parameters and the significance of the estimated errors further in Section 3.4.

Note that the resulting confidence intervals were determined after adding a systematic error of 0.025 mag in quadrature to the measurement errors listed in Tables 3 and 4. This accounts for measurement errors associated with the residual signal left by incompletely subtracted stars and galaxies along the target trajectory. Adding this offset affects only the confidence limits we find for the fitted parameters, not their best-fit values. The listed errors for H_g, H_r, and H_F606W do not include the calibration zero-point uncertainties (∼0.02 mag).

Figure 1 shows the minimum reduced χ² we find versus P, with the best fit occurring at P = 11.88190 hr with reduced χ² = 3.02 (239 observations and 7 fitted parameters). Note that there are many other χ² minima near the best fit (see right side of Figure 1). These are aliases with separation ΔP ∼ 0.027 hr, equivalent to one extra period in the long interval (438 days) between the two data sets most critical for fixing P (Gemini in 2017 July and DCT in 2018 October). The most significant secondary minimum occurs at P = 11.96254 hr with reduced χ² = 3.232. From Monte Carlo simulations, we find that these aliases yield a lower χ² than the true period less than 1% of the time. Hence, we reject this secondary minimum at the 99% confidence level.

Zoom In Zoom Out Reset image size

Figure 2 shows our best-fit sinusoidal light curve (solid line) as a function of rotational phase superimposed on the reduced observations. Here the observations are separated vertically by the date of observation, starting with the constructed HST values (M*, as described above), followed by the SOAR, Gemini, and DCT nights. See the figure legend for symbol explanations. Note that the bottom light curve is a second representation of the 2nd DCT epoch (2018 October 6) showing the assumed effect of increasing the brightness contribution from Thorondor (see discussion below). Also note that g-band data (red squares) are shifted relative to r-band data (black squares) by subtracting the fitted value for g–r listed in Table 5. All error bars include an added systematic error of 0.025 mag to account for residual, unsubtracted background flux along the target path in the images. The dashed lines show the nominally expected diminution from predicted mutual events at the relevant epochs, highlighted in green where the diminution exceeds 10% (discussed further below).

Zoom In Zoom Out Reset image size

It is clear by inspection of Figure 2 that our best-fit, single-harmonic sinusoid accounts for most of the signal variation we see from night to night (∼0.5 mag peak to peak), with the exception of the second epoch of DCT observations. We discuss the DCT observations further below. Excluding these DCT data, the most significant deviations from the fit occur on the first Gemini night (2017 July 22) over the phase range ω = 0.1–0.3. There the sinusoidal fit comes to gradual maximum whereas the observed light curve flattens and then dips by ∼0.2 mag (near ω = 0.25) relative to the fit. These deviations occur in both the g and r bands and are not an artifact due to unsubtracted background sources in the field along the target trajectory. The second peak (near ω = 0.7, covered only with the SOAR observations) does not appear to flatten significantly relative to the sinusoid. However, these data have larger scatter than the Gemini observations. There could be some flattening in the second peak that is unresolved.

With respect to the second-epoch DCT data, the best-fit sinusoid phases well with the observations, but the DCT data have a larger peak-to-peak variation (∼0.7 mag, compared to 0.5 mag for the model) and a fainter average signal (by 0.15 mag). The DCT light curve also flattens significantly from ω = 0.07 to 0.25 relative to the sinusoid. This mirrors a similar flattening in the Gemini light curve that occurs over the same ω range. On the assumption that the disparities in the DCT observations are the result of a diminished flux from Thorondor relative to earlier epochs (as might be expected from the ∼0.7 mag variability observed by HST), we computed the hypothetical light curve we would have measured if the flux from Thorondor were increased to its nominal value.

The resulting adjusted light curve appears at the bottom of Figure 2. Here we assume that the actual r-band magnitude of Thorondor at the time of the DCT observations was 25.4 mag (one magnitude fainter than the average F606W magnitude observed by HST). We brighten Thorondor by one magnitude to compute the adjusted light curve. It is apparent that the adjusted light curve is now in significantly better agreement with the best-fit sinusoid. The average brightness and the peak-to-peak variation of the two curves are nearly matched. Refitting the observations with this adjustment improves the reduced χ² from 3.02 to 2.39. Significant differences remain, however. Between ω = 0.0 and 0.2, the adjusted light curve rises more steeply (2.9 mag/period) than the sinusoid (2.0 mag/period). This phase range was not recorded at earlier epochs, so we cannot verify if this inconsistency was present at earlier epochs.

The second row of parameter values in Table 5 lists the new values for the light-curve parameters that result from fitting the adjusted observations. The changes are insignificant for P, A, H_F606W, and g–r. Given the improved χ², this validates our uncertainty estimates for these parameters. There is a small but significant change in β (from 0.00 to 0.05 mag/deg), which is coupled to small decreases in H_g and H_r (by 0.04 mag). These differences are a likely indication of our fitting uncertainty, which are larger than we have calculated from Monte Carlo estimates.

From the above comparisons we can conclude that the combined light curve for Manwë and Thorondor is periodic, with the peaks and valleys consistently occurring at the same rotational phase from the earliest HST observations (2007 July) to the latest ground-based observations (2018 October). For most epochs the light curve is sinusoidal to first order. The significant differences seen in the latest DCT data could be explained by Thorondor's long-term variability. In Section 3.6 we present a multicomponent model to better model the overall shape of the observed light curves.

Below we examine the light curves measured at each observed epoch to determine if a mutual event occurred. As will become evident, it is likely we detected none of the predicted events, or that our detections were at most marginal. The reasons could be that the time span of our individual light curves did not sufficiently cover the uncertainty window of the event times. More interestingly, it could be that the nominal event predictions, while based on the most likely assumptions, were unrealistic because of unusual shapes or sizes for Manwë and Thorondor.

Epoch 1 (SOAR data). As shown by the dashed lines in Figure 2, we would have detected a significant diminution in brightness of Manwë–Thorondor on the second SOAR night if the maximum diminution (0.66 mag), the event duration (∼12 hr or 1 full rotational phase), and the event midtime were as predicted. Even if the event midtime were early by 8 hr or late by 10 hr, the expected diminution would exceed 0.1 mag at the beginning or end of the observational window, respectively. However, given the close similarity of our light curves from both SOAR nights, it is clear that such a significant diminution did not occur. Taking into account only the uncertainty in the timing of the event, our nondetection had a likelihood of ∼40%.

Epoch 2 (Gemini data). Based on the above outcome, we planned our Gemini observations to detect a subsequent inferior event on the assumption that the actual midtime would be outside the ±1σ uncertainty window. Unfortunately, bad weather prevented us from scheduling our second night optimally. We therefore have no coverage in epoch 2 to verify a late event time. However, the first Gemini night would have overlapped the nominally predicted event (with the diminution exceeding 0.1 mag) if the event predictions were early by at least 10 hr.

Comparing the phased light curves from our two Gemini nights to each other and to the earlier SOAR nights, it again appears that no significant diminution occurred. As discussed in the previous section, there are some significant deviations on the first Gemini night relative to the best-fit sinusoid. Some of these deviations (e.g., the dip at ω = 0.25, or else the deeper-than-expected minimum near ω = 0.45) might result from a mutual event if both maximum diminution and duration were smaller than expected. We address these deviations again after presenting an improved rotational model for the system in Section 3.6. The main conclusion is that the mutual event did not occur earlier than expected if the depth and duration were as nominally predicted. This leaves open the possibility the event was late.

Epoch 3 (DCT 2017 November 9). The DCT observations scheduled 2017 November 9 would have detected a nominally predicted inferior event if it were early by 20 hr (2σ). Given the poor weather conditions and sparse observations, no conclusions can be drawn. The observations are marginally consistent with the best-fit sinusoid at the 1σ level.

Epoch 4 (DCT 2018 October 6). The DCT observations scheduled on 2018 October 6 would have covered a nominally predicted inferior event if it occurred earlier than predicted by 5–14 hr. However, the expected depth for this event was small (0.05 mag). As discussed above, the light curve we observed on this night deviates significantly from the best-fit sinusoid, with deeper minimum and a lower average brightness. If the event depth were much larger than predicted, some of these variances might be the outcome of the mutual event. To be consistent with the negative detections at earlier epochs discussed above, this would imply that the model used to predict the mutual events is grossly in error, with the events predicted at earlier epochs not occurring at all, or occurring with much lower diminution than expected. We return to this discussion later, after presenting an improved rotational model.

To better fit both the short- and long-term variations in the combined light curve of Manwë and Thorondor, we construct a model to predict the variability for a multicomponent system. Here we assume Manwë is a contact binary, rotating with period P_M and composed of two equal-sized lobes with matching, uniform albedos. This system is orbited by Thorondor, an independent body rotating with period P_T. The model determines the relative positions of the Sun, Earth, and all other bodies in the system and computes the total brightness variations caused by the mutual occultations and eclipses. In this paper, we assume a Lambertian scattering law for all surfaces.

As a starting point, we assume Manwë is a barbell-shaped contact binary composed of two equal-sized, elongated lobes connected by a narrow bridge (see Figure 3). As shown by Gnat & Sari (2010) and again by Descamps (2015), this is the shape expected for a fluid body with normalized rotational velocity, W = [3π/ρG]^1/2/P ∼ 0.30, where ρ is the density and G is the gravitational constant. As W approaches the critical value 0.2815, the narrow bridge connecting the two lobes vanishes and the hydrostatic solution becomes a pair of disconnected tidally locked bodies. Descamps uses shape models with W ranging from 0.30 to 0.32 to fit the light curves of the suspected contact binary (139775) 2001 QG298 (Sheppard & Jewitt 2004) and the main-belt contact binaries (624) Hektor and (216) Kleopatra. Based on this earlier work, we choose the shape model corresponding to W = 0.31. With the rotation period we measure for Manwë (P = 11.88 hr), this requires ρ = 0.80 g cm⁻³. This density value is near the nominal value measured for KBOs of comparable size (Stansberry et al. 2008), near the value (0.75 g cm⁻³) estimated by Grundy et al. (2014) for the combined system bulk density, and within the range of values estimated for other likely contact binaries (Thirouin & Sheppard 2017, 2018). Hence, this shape model is a good starting point for modeling Manwë's light curve.

Zoom In Zoom Out Reset image size

For simplicity and to speed our computations, we approximate each lobe of the barbell shape calculated by Descamps by a pair of contacting spheres, A and B, with diameters d_A and d_B, aligned with the long axis of the barbell. To set their relative scale, we refer to the shape parameters a and b (the half-length of the barbell and its maximum cross-sectional radius, respectively) given by Table 1 of Descamps (2015). Scaled to the radius of a sphere with the same volume as the barbell, a = 2.284 and b = 0.694. We set d_A to a value slightly less than the maximum cross-sectional diameter (1.8 × b) and d_B to span the gap between sphere A and the barbell's midpoint (a–1.8 × b). This yields d_B/d_A = 0.830. At maximum brightness viewed perpendicular to its spin axis, Manwë then has a surface area equivalent to a single sphere with diameter d_M = d_A [2(1+(d_B/d_A)²)]^1/2 = 1.838d_A and a volume equivalent to a sphere of diameter d_MVOL = d_A [2(1+(d_B/d_A)³)]^1/3 = 1.465d_A = 0.797d_M. With this shape, we are able to produce a light curve differing by less than 10% from the one generated by Descamps (2015) for the case W = 0.32, assuming the same scattering law.

To model Thorondor's contribution to the long-term variations of the light curve, we assume a spherical body with diameter d_T and a brightness that oscillates sinusoidally with period P_T (single-peaked) and amplitude A_T mag. This variation could be due to albedo variations across the surface (e.g., a dark surface feature rotating into and out of view). We further assume Thorondor's brightest visible hemisphere has an average albedo matching Manwë's albedo. This would be consistent with the matching colors observed for Manwë and Thorondor (Grundy et al. 2014) and for most KBO binaries (Benecchi et al. 2009). We recognize that Thorondor's shape is likely elongated and the dominant contribution to Thorondor's brightness variations. This is an important consideration because in mutual events with Manwë, Thorondor's shape is a significant factor affecting the short-term time dependence of the light curve. We therefore do not expect our simplified multicomponent model to accurately predict light curves during these particular events. At all other times, however, the shape we assume for Thorondor is less important. We assume a spherical shape only as a convenience to calculate Thorondor's peak brightness relative to Manwë. We revisit this assumption later. As long as the resulting variability is slow, periodic, and approximately sinusoidal, our model predictions will be valid. To avoid the complication of modeling the mutual events, we assume that no significant mutual event occurred during the observations, consistent with the conclusion of the previous section.

Despite our simplifying assumptions, the input parameters we need to completely specify our multicomponent model are numerous (see Table 6): diameters d_T and d_M, the heliocentric orbit of the tertiary system, the orbital elements for Thorondor's motion about Manwë's center of mass, Manwë's rotation period P_M and the orientation of Manwë's equator (ascending node Ω and inclination i referenced to the J2000 equatorial frame), and parameters P_T and A_T characterizing Thorondor's rotational light curve. We must also specify a reference epoch (ε_T) for Thorondor's rotational phase (chosen to be the date of the most recent minimum) and the rotational angle (μ) of Manwë's long axis relative to the ascending node of Manwë's equator at a given epoch.

Table 6.  Multicomponent Model Parameters

dM	188 km	Diameter of sphere with equivalent area as Manwë
dMVOL	150 km	Diameter of sphere with equivalent volume as Manwë
dT	108 km	Diameter of Thorondor
mM	1.41 × 1018 kg	Mass of Manwë
mT/mM	0.38	Upper limit to Thorondor/Manwë mass ratioa
pM	0.06	Albedo of Manwë
pT	0.09	Albedo of Thorondor
β	0.0 mag/deg	Linear dependence on solar phase angle
PM	11.88190 hr	Rotation period of Manwë
t1	2457625.8	Reference Julian date for absolute magnitudes
ρ	0.80 g m−3	Average density of the Manwë/Thorondor system
εT	2458398.0	Julian date of minimum brightness for Thorondor
Ω	1920	Ascending node of Manwë's equatorb
i	270	Inclination of Manwë's equatorb
μ	70	Rotation angle of Manwë's long axisc
$\langle {\alpha }_{{\rm{M}}}\rangle$	740	Average angle between Manwë's spin and Manwë-Earth direction
PT	309.03 days	Rotation period of Thorondor (single-peaked)
AT	0.55 mag	Light-curve amplitude of Thorondor (peak-to-peak)
Hg	7.48 mag	Absolute g-band magnitude of Manwë/Thorondord
Hr	6.66 mag	Absolute r-band magnitude of Manwë/Thorondord
HM	6.94 mag	Absolute F606W-band magnitude of Manwëd
HT	7.73 mag	Absolute F606W-band magnitude of Thorondord

Notes.

^aAssuming equal density, ρ, for Manwë and Thorondor. ^bJ2000 equatorial coordinates. ^cRelative to the ascending node of Manwë's equator at JD 2454300.0. ^dAt epoch t₁.

Download table as:  ASCIITypeset image

Fortunately, most of these parameters are known or well constrained. We refer to the JPL Horizons website⁶ for the heliocentric orbital elements. We refer to Grundy et al. (2014) for Thorondor's orbit about Manwë but we adjust the orbital period within its range of error (adding 0.020 day) so that no mutual events occur during any of the epochs covered by the observations. We fix P_M at the value we initially determine by fitting a sinusoid to the observed light curves (see Table 5). We set ε_T = JD 2458398.0, coinciding with the second DCT epoch, when we expect Thorondor to be near minimum brightness. We can also set d_T and d_M by adopting the nominal respective values, 92 and 160 km, estimated by Grundy et al. (2014). However, we scale these values up slightly by a factor 1.16 so that d_M = 188 km and d_T = 108 km. We then have an average density for the system, ρ = m_sys/[(π/6)([0.839d_M]³ + d_T³)] = 0.80 g cm⁻³, matching the density required by our shape model.

There are nine free parameters that remain to be determined by fitting our multicomponent model to the observations: the rotational orientation of Manwë (Ω, i, and μ), the light-curve parameters for Thorondor (P_T, A_T), the absolute magnitudes of the combined system in the observed pass bands (H_g, H_r), and the separate absolute magnitudes for Manwë and Thorondor in the F606W bandpass (H_M and H_T). We reference each absolute magnitude to the fixed epoch, t₁ = JD 2457625.8, chosen to be a time when Manwë and Thorondor are both near the peak of their separate light curves. We note here that Ω and i are coupled. All pairs of values that yield the same angle, α_M, between Manwë's rotational pole and the Manwë-Earth direction (or antidirection) at a given epoch will yield the same predicted magnitudes. Because the Manwë-Earth direction changes only by ∼10° during the total time spanned by the observations (including HST), α_M is nearly constant. The geometry provides little leverage for decoupling Ω and i.

To obtain the fit, we use the same χ² minimization procedure described above (see Section 3.3). Instead of using a sinusoid to described the model light curve, F(t), in Equation (1), we substitute the light curve computed using our multicomponent model. As before, we simultaneously fit all the observations listed in Tables 3 and 4, but we now separate the HST photometry for Manwë and Thorondor into two independent data sets. For each value of F(t), we therefore calculate offsets δF_M(t) and δF_T(t), which we add to F(t) to yield the predicted magnitudes for Manwë and Thorondor, respectively. We also normalize F(t) differently than before, defining F(t₁) = 0 mag at epoch, t₁. In our earlier sinusoidal fit, we normalize the model light curve so that F(t) = 0 mag at its average value. This leads to a relative offset for the computed values of H_g, H_r, and H_F606W. Also, because of this normalization, it is only the ratio of model parameters d_T and d_M that affects the χ² values of our fits. If our shape model is correct, however, then we know the density of the system. The sizes are then constrained by the known total mass.

Table 6 lists the resulting best-fit values for the free parameters and for the fixed parameters discussed above. The reduced χ² of the fit is 2.0 for 255 observations with 9 free parameters. The residuals to the fit are listed in the last column of Table 3. This is a significant reduction compared to our sinusoidal fit (χ² = 3.02). Figure 4 shows the best-fit model light curve compared to the observations. Except for the top curve, these are the same observations we present in Figure 2, again plotted versus rotational phase, ω(P,t) = (t–t₀)/P_M, with the value for t₀ listed in Table 6. The top curve shows the residuals to the fit of the magnitudes for Manwë and Thorondor separately measured with HST. Figure 5 shows the residuals for all the measured magnitudes plotted versus θ (averaged over bins of width Δθ = 001).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows that the new light-curve fit is in significantly better agreement with the observations than the sinusoidal fit shown by Figure 2. The model now matches the maximum of the light curve between rotational phases 0.1 and 0.4 (Gemini 2017 July 22 and DCT 2018 October 6) while simultaneously fitting the V-shaped minima appearing in the Gemini and DCT light curves. Taking into account the long-term rotational variability of Thorondor, the model is now able to match the relatively deeper minimum and overall decrease in brightness for the combined system observed by DCT on 2018 October 6. There are some observations of Thorondor near rotational phase ∼1.0 that have large scatter, and most fall well below the model prediction (magenta points in top curve). However, all these outliers correspond to the epochs when Manwë and Thorondor are separated by less than 008 (JD 2454399 and 2444766), making the deconvolution of the separate fluxes of Manwë and Thorondor unreliable. Because these observations have large errors, they have negligible influence on the model fit.

Note that we do not make an exhaustive search of parameter space to obtain this fit. We first sample the parameters on a coarse grid, and then search a finer grid only near the best-fit value of the coarse search. The particular orientation we obtain for Manwë's equator (Ω = 192°, i = 27°) corresponds to a pole orientation with α_M ranging from 690 to 789 over the time spanned by the observations in Tables 3 and 4. The pole's angle with respect to Thorondor's orbital pole is 125. However, as we discuss above there are other pole orientations that have consistent values for α_M and yield fits with only slightly higher χ². For example, we obtain χ² = 2.25 with Ω = 339° and i = 30°. The corresponding pole has a very different angle with respect to Thorondor's pole (555) while the α_M range (776–789) is not significantly changed.

Regarding the best-fit value for P_T, we find alternate solutions with P_T = 3.2359 and 331.13 days that yield nearly the same χ² as the best fit at P_T = 309.03 days. The shorter period is very close to an integer ratio (13/2) with Manwë's rotation period (0.495079 day), while the longer period is nearly an integer multiple (3) of Thorondor's orbital period (110.176 days). Both of these alternative values are likely aliases of the best-fit period generated by sampling bias. When fitting simulated observations that have the same observation dates and error bars as the real data but with magnitudes generated using our multicomponent model (assuming P_T = 309.03 days), we find these same aliases. In 94/100 simulated trials, however, the true period yields the best fit. Assuming our model is a realistic description of the system, we have 94% confidence that our fit has identified the true period. If the actual light curve for Thorondor is not sinusoidal, other period solutions are possible.

Table 5 lists the best-fit color, g–r = H_g–H_r = 0.81 ± 0.03 mag we obtain for Manwë–Thorondor assuming a sinusoidal fit. Our multicomponent fit yields nearly the same value (0.82). This is consistent with the value, g–r = 0.85 ± 0.06 mag, reported by Sheppard (2012). From inspection of Figure 2 or 4, we can also conclude that there is no significant deviation between the g- and r-band light curves. At the level of our g-band measurement errors (5%–10%), there is no dependence of color on rotational phase. Because Manwë's rotational variability is the dominant contribution to the light-curve variability on short timescales, we can conclude that Manwë's hemispherically averaged color is uniform with rotational longitude, although there might be a dependence on latitude. To the extent that color is a tracer of albedo (as it is on Pluto; Olkin et al. 2017), this also indicates a uniform hemispherically averaged albedo.

We derive an unusually flat slope to the solar phase curve for both Manwë and Thorondor from our sinusoidal fit (β = 0.0 ± 0.01 mag/deg) and we obtain the same value with our multicomponent fit. This flat slope is demonstrated by Figure 5, which shows that the residual magnitude to the multicomponent fit has no significant dependence on θ. While similarly flat phase curves for θ in the range 01–10 have been observed for icy satellites, the dwarf planets in the Kuiper Belt, and for the Haumea family members, most of the smaller and redder objects in the Kuiper Belt (with sizes and color similar to Manwë and Thorondor) have steeper phase curves (Rabinowitz et al. 2008; Verbiscer et al. 2019). The relatively flat phase curves for both Manwë and Thorondor require surfaces for both bodies with composition and/or scattering properties different from the general population (e.g., much higher albedo and larger ice fraction). Albedos much higher than we estimate from our multicomponent model (see Table 6) would require smaller sizes and higher densities than we have estimated for both bodies, inconsistent with the shape model we have assumed. At very small phase angles (θ < 01, near the lower θ range of our observations), the icy satellites exhibit strong opposition surges (Verbiscer et al. 2004; Buratti et al. 2011). Our observations do not extend to a sufficiently low phase angle to rule out this possibility for Manwë or Thorondor.

We note here that our analysis has not allowed for independent values of β for Manwë and Thorondor. It is possible that Thorondor has a steeper solar phase curve than the value we derive for Manwë and Thorondor together. The phase curve is sensitive to the rotational light curve we assume for Thorondor, which is not well constrained by the observations (see further discussion below).

We have shown that Manwë–Thorondor has a clear periodic component with a fixed period over the entire span of the observations. To first order (∼10% precision) and over a time span of a few days, the variability is sinusoidal. However, we obtain a better fit at all epochs (including fits to the separate Manwë and Thorondor photometry from HST) with a multicomponent model where Manwë is a barbell-shaped, strengthless contact binary rotating at the dominant period (with two peaks per rotation) while Thorondor varies sinusoidally with a single-peaked period of ∼309 days and amplitude ∼0.55 mag. It is remarkable that we obtain this fit by tuning only six parameters (Thorondor's rotation period and amplitude, Manwë's 2D spin orientation, and the rotational phases of Manwë and Thorondor). We must also make a small adjustment to the orbital period of Thorondor so that no eclipses occur, but no fine tuning is required for that purpose. It follows that our shape model is at least consistent with the observations.

If the shape model we have used is correct and Manwë is a strengthless body in hydrostatic equilibrium, then Manwë's density (ρ ∼0.8 g cm⁻³) is uniquely determined by its rotation period (P_M ∼11.88 hr). The equivalent spherical diameters for Manwë's volume and area (d_MVOL = 150 km and d_M = 188 km) and its mass (m_M = 1.41 × 10¹⁸ kg) follow directly from the known system mass (m_sys = 1.941 × 10¹⁸ kg) as long we have a good estimate for Thorondor's relative mass, m_T/m_M. As we discuss below, our model fit to the photometry suggests m_T/m_M = 0.38 (or less if Thorondor's shape is aspherical). It then follows that our estimates for d_MVOL and d_M are accurate to ∼15% precision. Given the absolute magnitude we measure at maximum brightness in the F606W band, H_M = 6.94, and assuming Manwë's V-band magnitude is fainter than the F606W-band magnitude by Δ_V,F606W = 0.383 mag (derived using synphot and Manwë's F606W–F814W color; see Benecchi et al. 2009), we can then derive an albedo for Manwë p_M = 0.06. This is within the nominal range of albedos known for KBOs the size of Manwë (Stansberry et al. 2008; Brown & Butler 2017).

Regarding the shape and rotation of Thorondor, we have no good constraints except that the brightness variability must be at least as large as we have measured (>∼0.55 mag) and the timescale of the variability must be ∼300 days if it is periodic. However, we do not have enough observations to verify periodicity. Our model fit does require that the surface-area equivalent spherical diameters of Manwë and Thorondor at maximum brightness have a ratio close to the values we assume, d_T/d_M = 108/188 = 0.57 (assuming matching albedos). If Thorondor is spherical and has the same density as Manwë, then we derive the relative mass, m_T/m_M = [d_T/d_MVOL]³ = 0.38. Here we use the relation d_MVOL/d_M = 0.797, which derives from the shape model we have assumed for Manwë. This estimate for m_T/m_M is an upper limit if Thorondor's shape is nonspherical because a sphere has the smallest surface area to volume ratio of any shape.

From Section 4.3, our assumption that Manwë is strengthless and in hydrostatic equilibrium allows us to determine d_M with better than 15% certainty. We can then accept our estimate for Thorondor's diameter, d_T = 108 km, as similarly accurate for a sphere with an albedo at maximum brightness matching Manwë's albedo. However, from the absolute magnitude we measure for Thorondor (H_T = 7.73 mag), and assuming the same F606W- to V-band transformation we use for Manwë (Δ_V,F606W = 0.383 mag), we derive a maximum-brightness albedo p_T = 0.09, significantly larger than the corresponding albedo we derive for Manwë (p_M = 0.06) This violates an implicit assumption of our model. It is not possible for both our model and our observations to be correct in this case.

The likely source of this conflict is our assumption that Thorondor's light curve is strictly periodic and sinusoidal over the entire time spanned by the observations. If we drop this assumption, we can artificially construct a light curve for Thorondor that yields a lower value for p_T. For example, if we alter our model to predict a brightness for Thorondor that is ∼50% higher during the time period spanned by the HST observations (2007–2013) compared to later observations (2016–2018), then H_T increases by 0.44 mag and p_T matches the value we derive for Manwë.

This brightness change would occur naturally if Thorondor is lenticular shaped (similar to the "Ultima" lobe of (486958) Arrokoth; Stern et al. 2019), rotating about its shortest axis (i.e., the lowest-energy rotation state) and if Thorondor's rotation pole is closely aligned with its orbital pole. The orbital pole's angle, α_T, with respect to Earth's line of site varies from 816 to 863 during the HST epochs, and from 901 to 914 during the later observed epochs. To calculate the resulting change in brightness, we assume Thorondor is a triaxial ellipsoid with semimajor axes, a > b > c. With a/c = 7.33, and b/c = 6.67, we have a flattened shape similar to the Ultima lobe, but with roughly half the relative thickness. Then from Equation (4) of Connelly & Ostro (1984), Thorondor's projected area, A, viewed from Earth is:

$$\begin{eqnarray*}&&\begin{array}{l}A=\pi {ab}\left[{\left(\sin {\alpha }_{T}\sin \omega /[a/c]\right)}^{2}+{\left(\cos {\alpha }_{T}\cos \omega /[b/c]\right)}^{2}\right.\\ \quad \,\,+{\left.{\left(\cos {\alpha }_{T}\right)}^{2}\right]}^{1/2}\end{array}\end{eqnarray*}$$

where ω is the rotational phase. This yields a mean visible surface area ∼40% higher during the HST epochs than the later observed epochs, nearly matching the long-term variability we require for Thorondor.

Although a lenticular shape would explain the long-term changes in Thorondor's brightness, it would not explain the ∼300 day variability we observe with amplitude ∼0.55 mag. However, such variability might arise from the precession of the orbital pole for this shape. From Equation (12) of Dobrovolskis (1995), the instantaneous rotational acceleration, d²ω/dt², for a triaxial ellipsoid on an eccentric orbit about a spherical primary at periapsis is:

$$\begin{eqnarray*}&&{d}^{2}\omega /{{dt}}^{2}=-K{(1-e)}^{6}-6\pi \gamma {T}^{-2}{(1-e)}^{3}\sin (2\nu )\end{eqnarray*}$$

where T is the satellite's orbital period, γ = (a²–b²)/(a² + b²), ν is the angle between the long axis and the direction to the primary, and K is ∼1.6 × 10⁻¹¹ yr⁻² for an icy body the size of Thorondor (scaled from the value estimated by Dobrovolskis for Nereid). The extrema for d²ω/dt² occur when ν is an odd multiple of π/4. With T = 110.176 days (see Table 1) and γ = 0.095 for the lenticular shape we have assumed for Thorondor, the extrema equate to ±0.0056 day⁻². Given Thorondor's slow rotation, this torque would be enough to alter Thorondor's rotational velocity by as much as 25% in one orbit. If Thorondor's rotational pole is tilted ∼125 from its orbital pole (see discussion above), then the maximum torque would precess the rotational pole by ∼70° about the orbital pole. Although the likely outcome over many orbits would be chaotic changes in Thorondor's rotation period and pole orientation, similar to the chaotic rotation initially suggested for Pluto's smaller moons (Showalter & Hamilton 2015), the resulting oscillations in α_T might yield the ∼300 day periodicity we observe for Thorondor at recent epochs.

Our multicomponent model does a good job of fitting the observed light curves on the assumption that no mutual events occurred during the observation windows. This is a plausible assumption given the large uncertainty in the event midtimes. However, we have explored alternative model fits for which significant mutual events do occur and significantly alter the model light curves. With fine tuning of the orbital parameters of Thorondor at values outside their ranges of uncertainty, it is possible to produce mutual event light curves that fit the observations approximately. If these solutions are correct, they would require significant secular changes to Thorondor's orbit over the time spanned by the observations. Additional HST imaging would likely be required to verify these orbital changes.

There is the additional possibility that mutual events occur during the observation windows but the resulting alteration of the light curves is marginal. For example, if Thorondor were to eclipse or occult the narrow connecting bridge between the lobes of Manwë, this would significantly reduce the expected depth and duration of the event. This scenario might be possible with an appropriate choice for the direction of Manwë's pole (while keeping the same angle with Earth's direction). If Thorondor has a lenticular shape that we view edge on during the mutual events (see discussion above), then this would further shorten the events and reduce their depths. It is conceivable that such an event is the cause of the dip observed in the light curve with Gemini on 2017 July 22 at ω = 0.25 (see discussion in Section 3.5). We plan to explore such scenarios in a subsequent paper.

Regarding the detectability of future mutual events, there is a possibility that the event season is either longer or shorter than nominally predicted. Grundy et al. (2014) assumed a spherical shape for Manwë. However, if the barbell shape we obtain for Manwë is correct, then the orientation of Manwë's rotational pole is a significant factor affecting both the duration of the events and the range of dates when the relative position of the Earth and Sun allow eclipses and occultations. Relative to the nominal diameter for the eclipse predictions (160 km) assumed by Grundy et al., Manwë may be larger by a factor of 2 in one orientation, but smaller by a factor of 2 in the perpendicular direction. The shape and orientation of Thorondor also affect the predictions. For favorable orientations, the event season may be a year or two longer than previously predicted.