key: cord-0601839-mjw0nz7j authors: Kobulnicky, Henry A.; Molnar, Lawrence A.; Cook, Evan M.; Henderson, Lauren E. title: A Bayesian Analysis of Physical Parameters for 783 Kepler Close Binaries: Extreme-Mass-Ratio Systems and a New Mass Ratio versus Period Lower Limit date: 2022-02-02 journal: nan DOI: nan sha: 1b51687d7f41acc3a254416fd7da60cf38175ad9 doc_id: 601839 cord_uid: mjw0nz7j Contact binary star systems represent the long-lived penultimate phase of binary evolution. Population statistics of their physical parameters inform understanding of binary evolutionary pathways and end products. We use light curves and new optical spectroscopy to conduct a pilot study of ten (near-)contact systems in the long-period (P>0.5 d) tail of close binaries in the Kepler field. We couple PHOEBE light curve models with Markov-Chain Monte Carlo analyses to compute Bayesian probabilities on inclinations, fillout factors, mass ratios, third-light fractions, and component temperature ratios. Mass ratios and third-light contributions measured from spectra agree well with those inferred from the light curves. Most binaries in the pilot study have extreme mass ratios q<0.32. At least eight are probable triples. A Bayesian analysis of all 783 Kepler 0.15 d0.5 d, as are systems with mass ratios near unity (q>0.8). There exists an empirical mass ratio lower limit $q_{min}$(P)~0.05-0.15 below which contact systems are absent, in agreement with a new set of theoretical predictions we obtain by modeling the evolution of contact systems under the constraints of mass and angular momentum conservation. Pre-merger systems should be found at at long periods and near this mass ratio lower limit, which rises from q=0.044 for P=0.74 d to q=0.15 at P=2.0 d. These findings support a scenario whereby mass transfer drives systems toward extreme $q$ and larger $P$ until the onset of the Darwin instability at $q_{min}$ precipitates a merger. The evolutionary paths of many stars culminate in a merger with a close stellar companion. As currently envisioned, these cataclysmic events occur after the stars undergo long geriatric episodes, exchanging mass as members of the ubiquitous population of contact binary systems. V1309 Sco became the prototype for stellar merger events when this 1.4 d contact binary exhibited an exponentially decreasing period, brightened, underwent a rapid evolution in light curve morphology, and erupted in a "luminous red nova" (Kulkarni et al. 2007 ) event similar to the 2002 V838 Mon eruption (Munari et al. 2002) , leaving only a single cool inflated star (Tylenda et al. 2011) . Stellar mergers may be as frequent as 0.2 yr −1 in the Milky Way (Kochanek et al. 2014; Howitt et al. 2020 ). The number of observed red novae in nearby galaxies is currently small, but the advent of wide-area deep sky surveys may precipitate detection of hundreds per year in the local universe. A comprehensive picture of the evolutionary sequence(s) yielding contact binaries and subsequent mergers is still elusive, but rapid advances are possible in this nascent era of all-sky time-domain datasets. Eggleton (2012) summarized a working theoretical sequence for the evolution of contact binaries terminating in a merger (see also Lucy 1976; Webbink 1976; Hilditch 1989; Stepien 1995; Yakut & Eggleton 2005) . The sequence commences with a wide (period P =months-1,000 yr) binary with main-sequence components orbited by a distant tertiary that induces Kozai-Lidov cycles (Kozai 1962; Lidov 1962) with tidal friction that shrinks the inner orbit to a few days, a point where magnetic braking can continue to rob the orbit of angular momentum. When the semi-major axis becomes small enough that the more massive primary star fills its Roche lobe, mass transfer to the secondary commences. Conservative mass transfer to the less massive component drives the orbit toward smaller separations and heats the secondary which fills its Roche lobe, eventually inducing mass transfer the other direction. A series of "thermal relaxation oscillations" drives the system deeper into contact and toward small mass ratios (q=M 2 /M 1 ). Once q reaches a critical threshold in the vicinity of 0.15, dictated by the point at which the tidally locked components' rotational angular momentum exceeds 1/3 their orbital momentum, the Darwin instability (Darwin 1893; Counselman 1973) instigates a rapid loss of angular momentum driven by tidal dissipation and non-conservative mass loss through the second Lagrange point. The runaway angular momentum loss culminates in a phase of common envelope evolution as the primary subsumes the secondary. Molnar et al. (2019) computed this critical ratio of angular momentum for a range of initial primary masses and system mass ratios, refining the mass ratio for the onset of the Darwin instability to a range of q=0.05-0.15 for systems with initial primary mass 0.6-2.0 M . This scenario predicts a paucity of systems with mass ratios more extreme than 0.10 and entails the requirement that contact binaries evolve from short-period (0.3-0.5 d) orbits toward ≈1 d orbits prior to coalescence. However, even if this sequence describes one dominant paradigm for contact binary evolution, alternate channels (e.g., dynamical interactions, formation driven by nuclear evolution) are likely to operate as well. Population studies of contact binaries lend some support to this "standard" scenario. W-type (Binnendijk 1970 ) contact binaries (those where the less massive secondary star appears to be hotter and produces the deeper minimum when eclipsed) are the most numerous and preferentially have shorter orbital periods in the 0.3-0.5 d range. By contrast, A-type contact binaries appear to have hotter primaries and tend to have longer orbital periods, lower mass ratios (Yakut & Eggleton 2005) and larger component radii, indicating a more evolved state (Mochnacki 1981) . Molnar et al. (2019) used a grid of MESA (Paxton et al. 2011 ) stellar evolution models to show that these observations are consistent with the expectations of mass and angular momentum conservation that drive contact binaries toward longer periods and low mass ratios. They also identified a small (N<100) class of long-period P ≈1 d contact binaries from among 22,400 close binaries in the 14-year OGLE photometric survey (Pietrukowicz et al. 2017 ) that exhibit large (|dP/dt| > 1.4 × 10 −8 ) negative period derivatives, consistent with being pre-merger candidates. Molnar et al. (2020) analyzed the light curves of 184,000 OGLE contact binaries to demonstrate an anti-correlation between mass ratio and period, supporting the evolution from W-type short-period toward A-type longer-period systems, a trend that we re-examine here using a sample of 783 short-period binaries having high-precision Kepler photometric data. High-quality photometry can provide excellent constraints on most system parameters in contact binaries and even in some detached binaries having large ellipsoidal light curve modulations. Fundamental parameters include periods, period derivatives, orbital inclination i, mass ratio q, fillout factor f , and ratio of stellar temperatures T 2 /T 1 . The fillout factor is a measure of the Roche lobe volume occupied and is defined in terms of the potential at the poles, Ω, and the critical potentials at the first and second Lagrange points, e.g., Prša (2018, eqn. 3 .67 and discussion therein). Fillout factors f <0 correspond to detached systems while f >0-1 correspond to contact systems. Figure 1 presents a graphical depiction of the Roche surfaces for fillout factors f=[0.05, 0.35, 0.65, 0.99] generated with PHOEBE for a q=0.25 system with a 1 M primary star. It illustrates the transition from nearly detached systems at f =0.05 (upper left panel) to the highly distorted secondary (nearly overflowing at the L 2 point) in the f =0.99 case (lower right panel). Figure 2 shows a grid of model contact binary light curves produced using PHOEBE for inclination angles i=[90 • , 60 • , 30 • , 10 • ] in rows top-to-bottom, fillout factors f =[0. 99, 0.70, 0.40, 0.20] in columns left to right, and mass ratios q= [0.95, 0.70, 0.35, 0. 10] coded by color as indicated in the legend. Here we adopt equal-temperature stars, T 1 =T 2 , since contact binaries have atmospheres in thermal contact, but actual systems can exhibit slightly unequal component temperatures (Hilditch et al. 1988; Yakut & Eggleton 2005) , sometimes as a consequence of spots (Barnes et al. 2004 ; Nelson et al. 2014 ). Phase φ=0 is defined, for purposes of this plot, when the more massive star is at superior conjunction. At large inclinations (top row, i=90 • ) the distinctive v-shaped or flat minima provide information on the inclination and mass ratio, while the overall amplitude of modulation, which decreases from left to right, provides constraints on f . Secondary minima are flatter and less deep as mass ratios becomes more extreme from q=0.95 toward q=0.10. At intermediate inclinations (second row, i=60 • ) the light curves become quasi-sinusoidal and more uniform, with the overall amplitude of modulation decreasing with decreasing q and f . At this inclination the more massive component still produces a deeper eclipse at φ=0. At still lower inclinations (third row, i=30 • ) the secondary minima (phase φ=0.5 when the less massive star is at superior conjunction) become deeper than the one at φ=0, more dramatically so as q becomes more extreme. The full amplitude of modulation is now 10%. At the lowest inclinations (bottom row) the minima at φ=0.5 become much deeper and wider than those at φ=0.0, most dramatically so at extreme q. The amplitude of modulation is also very small, on the order of 1%. One consequence of this reversal in timing of the deeper minimum with decreasing inclination is that the standard observational practice of locating the deeper minimum at φ=0 results in mass ratios M 2 /M 1 ≡q<1 being observed at i 50 • and q>1 for i 50 • -i.e., an eclipse of the more massive star always produces the deeper minimum (for equal temperatures). A fourth fundamental parameter, l3, the fraction of light from an unresolved third stellar component (either physically related or projected along the line of sight; not illustrated in Figure 2 ) serves to dilute the light curve modulation at all phases and alter the ratios of minima to maxima in a way that is partially degenerate with other parameters, including the temperature ratio, T 1 /T 2 . Figure 3 shows model light curves for contact binaries having i=60 • and f =0.5 for two mass ratios (q=[0.10, 0.25]) and three third-light fractions (l3=[0.0, 0.30, 0.90]), as coded by color and line style in the legend. Third-light contributions dilute the amplitude of modulation such that modest mass ratios with substantial third light become almost indistinguishable from more extreme mass ratios having small third light. For example, the q=0.10, l3=0.00 model (dashed black curve) is nearly (but not exactly) identical to the q=0.25, l3=0.30 model (solid cyan curve). Failing to diagnose third light correctly can lead to an erroneous mass ratio. The degeneracy is less severe at high inclinations and large fillout factors where light curve shapes are less ambiguous. With high-quality light curves, such as those afforded by Kepler, it is often possible to model and recover several or all of the principal system parameters. However, Figure 3 serves to illustrate danger in attempting to recover mass ratios solely from low-quality light curves when third-light contributions are unconstrained. If the third component has a spectral energy distribution substantially different from the contact binary, multi-color photometry can help break this degeneracy. However, in this work, we focus on the constraints afforded by high-quality single-band light curves. Additional free parameters such, as the ratio of stellar temperatures, T 2 /T 1 ( 1 for contact systems), and the possibility of inhomogeneities (e.g., hot or cool spots on the stellar surfaces) introduce additional signatures in the light curve that may be partially degenerate with other parameters but can nevertheless be modeled statistically, especially if multi-color light curves are available. Although detailed modeling of high-quality light curves can constrain system parameters of close binaries in many cases, phase-resolved kinematic measurements from spectroscopic data are required to measure individual component masses, total system masses, and constrain the structure and location of asymmetries in the systems such as hot or cool spots (e.g., techniques usually known as Doppler or Roche imaging; Vogt & Penrod 1983; Rutten & Dhillon 1994) . Spectroscopic data also serve to validate the results obtained from light curve modeling or provide mass ratios for detached binaries which can then be modeled to recover remaining parameters on the basis of light curves. Our goals in this contribution are 1) to critically assess the prospects for identifying extreme-mass-ratio contact binary systems from high-quality photometry in conjunction with state-of-the-art stellar binary model light curves, and 2) to test the Molnar et al. (2019) hypothesis that extreme-mass-ratio contact binaries are rare at periods exceeding ≈0.5 d and non-existent below a critical q threshold demarcating the onset of rapid stellar mergers. Section 2 describes the acquisition and reduction of new optical spectroscopic data of ten (near-) All models: f=0.5, i=60 Figure 3. PHOEBE model light curves for contact binaries with f =0.5, i=60 • , two different mass ratios (indicated by line style) and three different third-light fractions (indicated by color). Modest mass ratios with substantial third light closely resemble more extreme mass ratios with minimal third light. phases. The spectra directly provide the mass ratios and velocity amplitudes, informing light curve models that subsequently constrain the individual component masses and total system mass once the inclination is known. Furthermore, the spectral line velocity profile of the system serves to validate the remaining system parameters retrieved from light curve analysis alone (f , q, l3). Section 3 describes our application of the spectroscopic Broadening Function (Rucinski 1992) analysis techniques to determine the velocity profile of each of the ten contact systems in the pilot spectroscopic study. Section 3 also introduces our application of PHOEBE light curve models in conjunction with emcee (Foreman-Mackey et al. 2013) Markov-Chain Monte Carlo (MCMC) software to retrieve the Bayesian (posterior) probability distribution of system parameters. Section 4 employs Kepler spacecraft light curves in tandem with our spectroscopic data to measure the full set of system parameters for ten (near-)contact systems. We demonstrate the power of these combined datasets to vet state-of-the-art binary models while identifying extreme-mass-ratio systems possibly in the penultimate phase of evolution. We also use PHOEBE in conjunction with an MCMC analysis to sample the posterior probability distributions of system parameters as constrained by the data and obtain rigorous uncertainties on each. Section 5 extends use of these PHOEBE+MCMC retrieval tools to the entire set of >800 close Kepler binaries to obtain best-fitting and probabilistic system parameters in a uniform manner for an unbiased sample of unprecedented size and photometric precision. We investigate both contact and detached configurations for each system and identify the best-fitting configuration, where possible. Section 5 also investigates interesting statistical correlations among system parameters, providing insight regarding the evolutionary paths of close binaries. Section 6 provides a summary of substantial successes of the joint PHOEBE+MCMC approach and reviews the insights gleaned from the pilot study and the large-scale analysis of Kepler binaries that inform evolutionary scenarios for contact systems. Our intention throughout is to be pedagogical regarding some aspects of contact binary light curve analysis where we feel the modern literature is lacking and to be prescriptive in ways that help pave a path for large-scale analyses of binary systems in the age of vast time-domain datasets. We adopt the classical observational definition for eclipsing binaries that the primary star (mass M 1 , radius R 1 , effective temperature T 1 ) is the component that is eclipsed at the time of superior conjunction (t 0 ), producing the deeper eclipse. By this definition, the primary star is not necessarily the most massive or largest or hottest. This may lead to reported mass ratios, q=M 2 /M 1 , either greater than unity (corresponding to the W-type class of contact binaries in which the less massive component appears to be hotter) or q less than unity (corresponding to the A-type class of contact binaries in which the more massive component appears to be hotter). In some cases the primary minimum (orbital phase φ=0.0) and secondary minimum (orbital phase φ=0.50) have essentially the same depth, so the distinction between primary and secondary becomes ambiguous and arbitrary. During the analyses we are often interested only in the magnitude of the mass ratio without regard to which component is more massive. In such cases we invert mass ratios >1 to obtain what can be regarded as an absolute mass ratio, q a ≡[0-1]. 2. PHOTOMETRIC AND SPECTROSCOPIC DATA A small sample of ten long-period (P =0.55-1.38 d) contact binary stars were selected from the compilation of 2878 Kepler binary stars (Kirk et al. 2016 ) for spectroscopic observation. The targets were selected for having light curves consistent with contact or near-contact systems and longer-than-average periods. Figure 4 shows the distribution of period (P ) in days, effective temperature (T eff ), primary eclipse depth (p depth ), and light curve morphology parameter (morph) 1 for all Kirk et al. (2016) binaries with periods between 0.2 d and 4 d (black points) and our targets (red star symbols). Kirk et al. (2016) note that morphology parameters between 0.5-0.7 correspond to semi-detached systems while 0.70 and higher belong to contact systems and ellipsoidal variables (i.e., tidally deformed detached systems). The dense locus of points in the lower left panel forming a linear trend of increasing temperature with period marks the population of close binaries with main-sequence components. The correlation reflects the relation between temperature and radius for main-sequence stars. This trend becomes less pronounced above about 7000 K, reflecting the small number of high-mass (M 2 M ) hot stars observed by Kepler, a consequence of both the stellar initial mass function and possible selection biases imposed by the Kepler Input Catalog (Brown et al. 2011; Batalha et al. 2010, KIC) . These hotter stars preferentially have morphology parameters greater than about 0.75 (appropriate to contact and ellipsoidal variables) and primary eclipse depths reflecting the full range of values, as expected from a random distribution of inclination angles. At periods longer than about 0.5 d the distribution of morphologies for Kepler binaries bifurcates into a lower branch reflecting detached morphologies and an upper branch indicating contact and ellipsoidal variables. Another narrow locus of points in the upper left panel near p depth 0 stretching from short toward long periods reflects the population of statistically numerous grazing-eclipse systems. Our target sample (red star symbols) has morphology parameters ranging from 0.74 to 0.95, relatively large primary eclipse depths, and periods ≈ 1 day, placing them on the long-period tail of the distribution. Table 1 lists the ten objects in the pilot study by Kepler Input Catalog identifier (column 1), mean Kepler band magnitude (column 2), stellar effective temperature from the KIC (column 3), light curve morphology parameter (column 4), orbital period (column 5), and reference time (t 0 ) of superior conjunction from Kirk et al. (2016) with updates in this work (column 6). We assembled calibrated Kepler photometry on the targets available from the the public MAST 2 archive as cleaned and detrended by Kirk et al. (2016) 3 . Data were generally available from a majority of the Kepler operational quarters, from as few as seven to as many as 17, yielding tens of thousands of measurements in the broad Kepler bandpass spanning four years, from 2009 May through 2013 May. We determined mean periods for each system using a Lomb-Scargle periodigram analysis 4 which were found to be in good agreement with those tabulated by Kirk et al. (2016) . These periods were then used to fold the light curve and calculate a time of superior conjunction, t 0 -the reference time centered on the deeper minimum in the light curve. In the manner of Molnar et al. (2017) we fit the full light curve using a sum of three-six Fourier components to define an analytic function characterizing the mean light curve. This function was then used to fit subsets of the light curve, one Kepler quarter at a time, to measure time-dependent phase shifts that indicate variations in the times of minima or maxima. Such variations, termed "eclipse timing variations", may indicate a changing orbital period or a light travel time delay that results from an orbit about a third body. In a . Distribution of orbital period, effective temperature, light curve shape morphology parameter, and primary eclipse depth for binary stars in the Kepler field (Kirk et al. 2016 , black points) and our ten targets (red stars). few cases where suitable data were available and helpful, we added recent 2018-2019 photometric measurements from the Zwicky Transient Factory (ZTF, Smith et al. 2014 ) as a way of extending the time baseline. We obtained optical spectra on the targets near each of the two quadrature orbital phases (φ=[0.25,0.75]) with the longslit spectrographs at the Wyoming Infrared Observatory (WIRO) 2.3 meter telescope and/or the Apache Point Observatory (APO) 3.5 meter telescope. At APO we used the Double Imaging Spectrograph (DIS) red arm with a 1200 line mm −1 grating in first order to acquire spectra over the range 5800-6900Å. The spectral resolution was ≈3000 using slit widths of 0. 9 or 1. 2 at a reciprocal dispersion of 0.62Å pix −1 . The wavelength calibration was performed with the aid of HeNeAr lamp spectra obtained close in time to the target exposures and has an RMS of 0.06Å, or about 3 km s −1 at 6500Å. Spectra acquired at WIRO used the Longslit spectrograph with a 1. 2 slit and a 2000 line mm −1 grating in first order to cover 5400-6700Å at a reciprocal dispersion of 0.61Å pix −1 and resolution R 4000. The CuAr lamp exposures acquired after each science exposure provided wavelength calibration to an RMS of 0.019Å (∼1 km s −1 ). At both observatories, target exposures times varied from two×600 s to four×600, depending on seeing and source brightness, yielding spectra with continuum signal-to-noise ratios (SNRs) of 40:1-100:1 near 6500Å in the final combined spectra. Averaging spectra over ≤40 minutes introduces a small amount of phase smearing which is small compared to the rotational broadening of the tidally-locked short-period stellar components. Spectra were reduced using standard techniques in IRAF (Tody 1986) , including 1D spectral extraction and local sky background subtraction, flat fielding with quartz continuum lamps, and transformation to the Heliocentric velocity frame of reference. Observations of radial velocity standard stars confirm that the velocity calibration is precise to ±6 km s −1 between observatories and epochs, with deviations primarily attributed to variable placement of a target within the slit-an inevitable limitation of slit spectrographs. We analyzed the optical spectra using a custom python version of the broadening function algorithm (BF) described by Rucinski (1992 Rucinski ( , 2002 to recover the velocity profile of each contact binary at each quadrature phase. The BF code performs a true linear deconvolution of a broadened stellar spectrum given a narrow-lined spectral template of the appropriate effective temperature. It is superior to cross-correlation methods when used to separate the blended components of close binary systems having similar temperature (Rucinski 1999) . The resulting BF for a contact binary is the light-weighted velocity profile of the combined system at the time of observation, as broadened by the instrumental profile (which is 75 km s −1 FWHM, smaller than the ≈200 km s −1 rotational profiles) and any temporal broadening from finite exposure durations. The BF can also reveal the signatures of third components in the spectra of rotationally broadened binary systems (e.g., D'Angelo et al. 2006) . Our BF analysis used as a template a high-resolution high-SNR spectrum of the appropriate T eff (we adopt solar metallicity and log g=3 models for inflated stars; the exact choice in inconsequential for our purposes) from the PHOENIX model atmospheres (Husser et al. 2013) . However, gross mismatches in T eff between the spectra and the template of more than about 1000 K produce negative "bowls" on both sides of the BF peak. Good matches between the data and the template produce the largest BF amplitudes, serving as a check on the suitability of the effective temperature listed in the Kepler Input Catalog. We compared BFs resulting from the full spectral coverage (usually ≈5450-6650Å) and from a spectral subregion that excludes Hα, the strongest single spectral feature. We found that BFs are consistent with each other regardless of spectral regime, but they have larger SNRs when Hα is included. In principle, low levels of Hα emission in the core of the line could lead to a skewed BF because of the large weight of this feature in the spectrum, but we found no evidence for such effects. For the ten systems studied spectroscopically (Section 4) and the entire ensemble of nearly 800 candidate contact binaries from the compilation of Kirk et al. (2016) (Section 5) we modeled the light curves and velocity curves using the binary modeling code PHOEBE 2.2 5 in conjunction with the Markov-Chain Monte Carlo code emcee (Foreman-Mackey et al. 2013) to explore the posterior probability distributions of system parameters. We adopted the period (P ) and time of superior conjunction (t 0 ) from Kirk et al. (2016) , except in a few cases where we recomputed slightly different values using a subset of the Kepler data. We adopted T 1 from the Kepler input Catalog (Brown et al. 2011) , using 6200 K if no value was listed. The models are insensitive to the adopted T 1 because the limb darkening coefficients vary only modestly over the 4500 K< T eff <7500 K range of contact binaries of concern here (e.g., limb darkening tables of van Hamme 1993; Claret 2004) . Given the short periods of these systems, we set zero orbital eccentricity (e=0) as a fixed parameter. We used the mean Kepler passband in the model light curve and the default Castelli & Kurucz (2003) stellar atmosphere models with limb darkening coefficients interpolated from these models as implemented in PHOEBE 2.2. All systems were modeled using four different computational approaches. • Fixed-temperature-ratio (T 1 =T 2 ) model-a contact binary geometry with equal temperature components (T 2 /T 1 =1) and four free parameters: inclination i, fillout factor f , mass ratio q, and third-light fraction l3. We ultimately concluded that this model is too restrictive, leading to incorrect solutions as the fitting process is forced to alter other model parameters to compensate for lack of flexibility in T 2 /T 1 . • Variable-temperature-ratio model-a contact binary geometry with five free parameters: inclination i, fillout factor f , mass ratio q, third-light fraction l3, and temperature ratio 0.7 M 1 . The model velocity semi-amplitude of the more massive and larger contact binary component is 4.2 km s −1 , consistent with the single peak at nearly constant velocity in the BF. The extreme mass ratio and inclination implies that the less massive component is faint and would be hard to detect (indeed, it is not detected) in the spectral profile despite its large velocity amplitude, especially if the third-light contribution to the BF is substantial. 15 The green x's in Figure 7 show the theoretical line profile generated by the best-fitting PHOEBE model at the observed orbital phases after convolution with the instrumental spectral profile. 16 The normalization of both the BF and the theoretical line profile is arbitrary, so they have been scaled to approximate the suggested l3≈0.7 third-light contribution which dominates the BF. It is not possible to draw conclusions about the luminosity of a third body on the basis of the BF since it is expected to be blended with the profile of the contact binary. In this system, only one component is apparent in the BF. We are unable to distinguish whether this is the brighter component of the contact binary or a third component. 17 The small offset of about 10 km s −1 between the theoretical profile and the BF at . Broadening function of KIC04853067 (black dots) showing a single peak with a small velocity difference between the two epochs. Red curves depict two Gaussian components fit to the data, while the blue curve is the sum of those components. Green x's show the theoretical line profile generated by PHOEBE for the best-fitting contact binary system parameters as convolved with the instrumental spectral profile. The normalization of both curves is arbitrary. The stellar components are viewed at low i and are unresolved. The BF may be dominated by a third component that is not represented in the model line profile function. This system may also be a detached configuration. φ=0.18 is consistent with epoch-to-epoch wavelength calibration uncertainties, but it could also include a contribution from orbital motion of the contact binary about a putative third star implied both by the O − C analysis and the ensuing third-light analysis. The inclination of the third body's orbit is unconstrained. To investigate the power of the light curve to constrain the system parameters in the absence of kinematic measurements, we performed a PHOEBE+MCMC analysis for each of the three competing models. Figure 9 shows contour plots (2D) and histograms (1D) of the relative probabilities of the four free parameters used in the fixedtemperature-ratio contact binary models (left panel) and the five free parameters of the variable-temperature-ratio models (right panel). The inner and outer contours enclose 39% and 86% of the samples, respectively, corresponding to 1σ and 2σ levels for a two-dimensional normal distribution. The simulations show that all parameters are unimodal and well-constrained. The one-dimensional histograms are notably non-Gaussian in all parameters. The most probable parameters and 1σ one-dimensional uncertainties of the fixed-temperature-ratio models-as represented by the 16th/50th/84th percentiles on each parameter-are cos i=0.912/0.928/0.943 (i 18 • ), f =0.15/0.36/0.69, log q=0.733/0.928/0.949, l3=0.75/0.79/0.83. The parameters of the best-fitting variable-temperature-ratio model are consistent with these ranges, with the exception of fillout factor which is much larger than the fixed-temperature-ratio models. Visualizations of 100 parameter sets randomly selected from the MCMC ensemble closely follow the bestfitting model depicted in Figure 8 , providing assurance that the walkers sampled the allowed parameter space in a suitably probabilistic manner. The best-fitting and the most probable solutions require low inclinations, highly uncertain fillout factors, extreme mass ratios, and significant third-light contributions, consistent with the eclipse timing variations depicted in Figure 6 . Figure 10 presents two-dimensional contour and one-dimensional histogram distributions for each of the six free parameters resulting from the detached model MCMC analyses. The permitted parameter range is large for many of the parameters except inclination. Nevertheless, KIC04853067's light curve provides meaningful constraints on the majority of key system parameters. The RMS of the best-fitting solutions is nearly identical to that in Figure 8 for the Figure 8 . Light curve and velocity curve of KIC04853067, along with best-fitting contact binary models (magenta: variabletemperature-ratio; green: fixed-temperature-ratio). Individual points represent Kepler photometric data and residuals (upper and middle panels) or spectroscopic measurements (lower panel). Labels denote 16th (best)50th 84th percentile values for the principal model parameters. The components' line profiles are likely blended and contaminated by third light so that individual velocity amplitudes cannot be measured for this low-inclination extreme-mass-ratio binary. contact configurations. The most probable T 2 /T 1 here is near unity, consistent with a contact configuration wherein the stellar components share an atmosphere. The most probable mass ratio near q=10 (q a =0.1) is extreme and both components are close to overflowing with R 1 /R 1max >0.9 and R 2 /R 2max ≈0.9. In such configurations, contact and detached systems become indistinguishable as the larger overflowing or nearly overflowing component dominates the light curve modulation. In conclusion, KIC04853067 is a P ≈1.3 d extreme-q candidate system exhibiting a well-measured low-amplitude light curve that can be modeled equally well as a contact or detached configuration. This is consistent with its large light curve morphology parameter of 0.88, approaching the regime of ellipsoidal variables. Either configuration requires a low inclination angle and both stars near overflowing. The system's velocity profile is kinematically unresolved, consistent with low i and large third light. Systematic O − C variations and a preference for significant l3 in the Bayesian analyses of the light curve indicate that KIC04853067 is a triple system in both the contact and detached models. The case of KIC0485306 is a cautionary tale of how even high-precision single-band light curves can fail to discriminate between contact and detached geometries when mass ratios are extreme or inclinations are small. . Two dimensional contours and one-dimensional histograms depicting the posterior probability distribution for combinations of free parameters for the fixed-temperature-ratio models (left) and variable-temperature-ratio models (right) used to model the KIC04853067 light curve. The inner and outer contours enclose 39% and 86% of the samples, respectively, corresponding to 1σ and 2σ uncertainties for normal distributions, but note the significantly non-Gaussian histograms in each parameter. Degeneracies between temperature ratio and other parameters are prominently visible by the elongated contours in the right panel. . Two dimensional contours and one-dimensional histograms depicting the posterior probability distribution for combinations of six free parameters (plus two ancillary derived parameters R1/R1max and R2/R2max) for the detached configuration models of KIC04853067. Model parameters remain well-constrained and most probable values overlap with the fixed-temperature-ratio models indicating low inclination, extreme q, and significant third light. Figure 11 shows the folded Kepler light curve of the P =0.99 d system KIC04999357. The Fourier components used to approximate the mean light curve (red curve) are given in the Figure. Primary and secondary eclipse depths are very nearly equal, as are the maxima between eclipses. The continuous modulation is that of a classical contact binary. KIC04999357 exhibits no evidence for O − C variations over the duration of the Kepler mission, with an RMS variation of <0.2 min. Accordingly, we adopt a linear ephemeris. Figure 11 . Folded Kepler light curve of KIC04999357 and mean light curve shape (red curve) generated from the sum of low-order Fourier coefficients. Figure 12 shows the BF (black dots) for two epochs of spectroscopy at orbital phases φ=0.24 and φ=0.77. There are two clear components. We fit the BF using a two-component Gaussian 18 function (red curves) which provide an estimate of the component velocities at each quadrature phase, and thereby, a mass ratio and a systemic velocity. The sum of the two Gaussian components (blue curve) provides a good match to the BF data. The ratio of areas of the two Gaussian components implies a luminosity ratio L 2 /L 1 of about 0.28 (1:3.5) and a radius ratio L 2 /L 1 = R 2 /R 1 =0.56. The radial velocities of the components yield velocity semi-amplitudes of K 1 =41 km s −1 and K 2 =187 km s −1 , implying a mass ratio near q=0.22 and a systemic velocity of γ = −60 km s −1 . However the velocity of the secondary at φ=0.24 is poorly measured, as it appears considerably fainter and less defined than at φ=0.77. A second spectrum obtained at φ=0.27 on 20170711 shows the same deficit near the expected peak of the secondary's BF. This deficit is present regardless of spectral range used in the BF analysis, effectively ruling out the possibility of emission in any one spectral feature or uncorrected Telluric absorption affecting the spectra. We find that introducing a cool spot on the trailing face of the secondary star in the PHOEBE models can roughly reproduce this deficit in the BF. The deficit here, and in some subsequent examples, bears some similarities to that of the secondary in extreme-q system AW UMa which Rucinski (2015) interprets as indicating the presence of an accretion disk or flow of matter to/from the secondary. Figure 12 . Broadening Function of KIC04999357 showing two velocity components. The black dotted curve is the BF, the red curves are two Gaussian components fitted to the BF, the blue curve is the sum of the Gaussian components, and the green x's represent the model line profile of the best-fitting PHOEBE model. The normalization of both the data and the model is arbitrary. KIC04999357 is the only one among our pilot sample to have large astrometric uncertainties, as judged by the Renormalized Unit Weight Error 19 (RUWE=3.00) in the Gaia DR2 dataset (Gaia Collaboration et al. 2016) . RUWE is a metric that designates a poor single-star astrometric solution when RUWE>1.4 and is often used to flag candidate multiple-star systems (e.g., Kervella et al. 2019) . Contact binaries having semi-major axes of several solar radii would not yield poor Gaia astrometric solutions on their own, but tertiaries orbiting contact binaries at distances few AU could produce large RUWE values. Hence, there is evidence for a third component in this system or along the line of sight. Figure 13 shows the folded photometric data (upper panel) and the spectroscopic radial velocity data (lower panel), with the best-fitting variable-temperature-ratio PHOEBE model (without irradiation effects) overplotted using solid magenta curves. Text within the panels also gives the 16th/(best)50th/84th percentile ranges from MCMC simulations. The best-fitting model (RMS=0.0018) requires i=67.9 • , f =0.70, q=0.17 (in reasonable agreement with the q=0.22 from the broadening function), l3=0.35, and T 2 /T 1 =1.02. This model is superior to the best-fitting detached models which all have larger RMS and require R 1 /R 1max >0.99 and R 2 /R 2max >0.99, indicating a contact configuration. The contact model yields a ratio of component radii R 2 /R 1 =0.48 (close to the value inferred above from the BF), R 1 =2.68 R , and R 2 =1.29 R . Adopting this inclination allows the computation of component masses M 1 =1.24 M and M 2 =0.21 M . The line profile function produced by this model (green x's in Figure 12 ) shows a reasonable match to the BF when normalized to the fainter component. However, the amplitude of the brighter component in the model is smaller than the peak in BF. The difference between the data (black dotted curve) and model (green x's) can readily be explained by the light of a third component amounting to ≈10% of the total system light, assuming the third component has a similar temperature to the contact binary. As the present release of PHOEBE does not support a full inclusion of tertiary components in the line profile, we cannot model the contribution of the putative third star to the BF here or in subsequent examples where third components are probable. The presence of a third star is likely to skew the primary peak in the BF, making measurement of the primary star's velocity uncertain. The middle panel of Figure 13 plots the model residuals (model minus data) of the best fitting variable-temperatureratio light curve (magenta curve) and the best-fitting light curve when irradiation and reflection effects are included in the model (green dashed curve). The best-fitting model with irradiation has a modestly smaller RMS (0.0014 versus 0.0018), but the 16th/50th/84th percentile ranges are very similar to those listed in the top panel and the best-fitting parameters are all nearly identical. Residuals are smaller, primarily in the region centered on the primary eclipse (φ=0), as expected based on the comparisons presented in Horvat et al. (2019) . While an attempt at proper treatment of irradiation effects improves the best-fitting models in this well-measured system having small intrinsic dispersion in the data, the impact on the Bayesian constraints on each parameter is minimal. Figure 13 . Folded light curve and velocity measurements of KIC04999357, along with best-fitting PHOEBE model light curves, residuals (model minus data), and velocity curves without irradiation effects (magenta curves) and with irradiation (green curves). Figure 14 shows the results of the MCMC analysis for KIC04999357 stemming from the nominal T 2 ≈T 1 models. The contours and peak density of points indicates that the most 16th/50th/84th percentile system parameters are cos i=0.34/0.38/0.43, f =0.52/0.67/0.81, log q=−0.90/−0.79/−0.67 (q=0.18, in general agreement with the BF), l3=0.16/0.31/0.43, and T 2 /T 1 =0.99/1.00/1.01, consistent with the best-fitting parameters given above. All five of the model parameters are well-constrained. The Monte Carlo simulations allow for a modest l3≈0.31, not inconsistent with the contribution inferred from the BF discussed in connection with Figure 12 previously. In summary, KIC04999357 (morph=0.84) is a possible triple system where the inner contact binary contributes most of the luminosity and has a small mass ratio q=0.17-0.22. PHOEBE models with and without irradiation effects produce convincing fits to the data with a component temperature ratio near unity. KIC04999357 may be an evolved binary that has exchanged mass, on its way toward the Darwin instability limit. The evidence for a third componentfrom the large GDR2 RUWE values and the excess in the BF near the systemic velocity-is consistent with the idea that Kozai-Lidov cycles initially play a role in bringing the inner components into contact. The putative third body must lie at a large separation from the contact binary in order that it produce a large astrometric RUWE and not produce detectable O − C variations. Figure 15 displays the folded Kepler light curve of the P =1.08 d system KIC06844489. The mean light curve displays a 10% semi-amplitude, very little dispersion, and is well characterized by the Fourier components labeled in the Figure. Only seven quarters of Kepler data are available. The secondary minimum is less deep than the primary minimum. . The red solid curve shows a best-fit parabola, representing a quadratic ephemeris with a period derivative of dP/dT = −4.66×10 −9 -a constantly decreasing period. This model is a poor match to the O − C data. The dashed green curve shows a sinusoidal model with amplitude A=12.1 min, P =11.7 yr, and t 0 =2455691 (the peak of the sine curve-the time when the contact binary is at superior conjunction relative to the foreground third component). This model, predicated on the hypothesis of a third body in orbit with the contact binary, provides an excellent match to the data. The semi-amplitude of this curve gives the light crossing time across the projected semi-major axis of the contact binary about the system Barycenter, t=asin(i outer )/c. The implied minimum semi-major axis of the contact binary's orbit is 1.45 AU. A hypothetical 0.34/sin(i outer ) M tertiary in an 11.7 yr orbit can explain the O − C data. The broadening function of KIC06844489 in Figure 17 shows two components, one much smaller than the other. The spectroscopic data at phases φ=0.25 and φ=0.76 (assuming a linear ephemeris) were obtained 1.5 months apart. The primary BF component has a distinctly non-Gaussian line profile, as if it may be a blend of unresolved subcomponents. The ratio of component velocities in a two-component fit indicates an extreme mass ratio near q=0.10. However the presence of a bright third component spectrally blended with the primary would result in the K 1 amplitude and q being underestimated. Given the evidence for a third body from the large O − C variations, we fit the BF using three components, fixing the velocity of the tertiary near −25 km s −1 . The sum of the three components (blue curve) shown by the three red curves in Figure 17 provides a good overall fit to the BF, with the third component contributing 60% of the total light (under the assumption of a similar temperature to the contact binary). The implied mass ratio is q 0.25, but the velocity of the fainter component is rather uncertain. There is considerable degeneracy between the strengths of the primary and tertiary, rendering the velocity of the primary highly uncertain. As in KIC04999357, the secondary's peak near φ=0.25 is smaller than at the other quadrature phase. Figure 18 shows a light curve and velocity curve of a best-fitting variable-temperature-ratio PHOEBE model. A model with i=73.0 • , f =0.98, q=0.31 (in good agreement with the BF), l3=0.63 (consistent with the deficit in the primary's BF peak in Figure 17 ), and T 2 /T 1 =0.98 produces an excellent fit (RMS=0.0010). Models implementing irradiation effects require nearly identical system parameters and do not improve the fit. The residuals in the middle panel are small and show no systematic offsets near phases 0 or 0.5 where irradiation effects would be expected to create the largest signatures. At this inclination the measured component velocities imply masses M 1 =1.77 M and M 2 =0.51 M . These masses dictate component equivalent radii R 2 /R 1 =0.64, R 1 =3.19 R , and R 2 =2.06 R . The best contact models fit the light curve well but produces a line profile (green x's in Figure 17 ) having primary star amplitude smaller than the primary BF peak, suggesting a need for a third component. By contrast, the best-fitting detached model has a much larger RMS and indicates the primary is nearly overflowing with R 1 /R 1max ≥0.99. Accordingly, we consider this to be a bona fide contact system with nearly equal temperature components and a large fillout factor. Monte Carlo simulations of the KIC06844489 system using the T 2 ≈T 1 model yield the probability distributions in Figure 19 . The most probable values are all well-constrained at cos i=0.292/0.311/0.338 (i 73 • ), f =0.79/0.91/0.97, log q=−0.58/−0.52/−0.46, l3=0.57/0.60/0.62, and T 2 /T 1 =0.96/0.97/0.98. The Figure shows some degeneracy between cos i and f in the sense that larger cos i (smaller inclination) requires smaller fillout factors. There is also degeneracy between f and T 2 /T 1 , where larger f demands larger T 2 /T 1 . Evidence for a substantial third-light component is supported by the MCMC analysis, consistent with the O − C residuals in Figure 16 and the excess in the BF near the systemic velocity. The minimum (i=90 • ) third-body mass implied by the O − C analysis of M 3 =0.33 M is not capable of producing 63% of the system light if it is a main-sequence star. Orbital inclinations of the M 3 -(M 1 +M 2 ) system of i outer <20 • would result in M 3 1.5 M , allowing a main-sequence F star tertiary to contribute substantially to the system light in the amount suggested by the Monte Carlo simulations. In summary, the evidence indicates KIC06844489 is a contact system (consistent with its morph=0.84) with a luminous third component that creates large O − C systematics well-modeled by a sine function. The most probable mass ratio of the contact binary system is small at q 0.30. It is a somewhat massive (M tot =2.28 M ) long-period contact binary system that may be nearing the putative Darwin Figure 21 shows the broadening function of KIC08913061 at phases φ=0.36 and φ=0.71. There are two distinct peaks in the BF with the small peak at positive velocities at φ=0.36, indicating a mass ratio q<1. The component velocities yield K 1 =36 km s −1 , K 2 =202 km s −1 , q=0.18, and systemic velocity γ=−72 km s −1 . As in the previous two examples, the secondary's BF peak is again smaller at φ=0.36 than at the opposite quadrature phase. We performed a PHOEBE Monte Carlo analysis of the KIC08913061 system, finding that the most probable system parameters involve large inclinations i>70 • , rather extreme mass ratios, and T 2 /T 1 1. The best-fitting variabletemperature-ratio model shown in Figure 22 (upper panel, magenta curve) provides a good fit (RMS=0.0042) to the light curve with i=81.0 • , f =0.50, q=4.2 (secondary more massive than primary), l3=0.39, and T 2 /T 1 =0.86 regardless of whether irradiation effects are included. Such a large temperature difference between the stars seems inconsistent with a contact system. Furthermore, this mass ratio is grossly inconsistent with the kinematic data which require q<1! As an alternative we considered a fixed-temperature-ratio model-T 2 /T 1 =1. The best solution (red curve) is i=90 • , f =0.99, q=0.28 (roughly consistent with the kinematic data), and l3=0.44, but this model produces a less satisfactory fit to the light curve (RMS=0.0141, about 3.5 times larger than the variable-temperature-ratio model), with large residuals (middle panel) near phases 0.0 and 0.5. Enabling irradiation effects does not improve model. Detached models fit the data less well than the variable-temperature-ratio contact models and require both components to have R 1 /R 1max >1 and R 2 /R 2max >1, indicating that a contact configuration is preferred. As a third scenario we tried a T 2 /T 1 =1 contact model with a hot spot on the primary star. 20 This model (green curve) produces a superior fit to the light curve (RMS=0.0029) while also matching the kinematic data for parameters i=77.0 • , f =0.46, q=0.16, l3=0.025, T spot /T 1 =1.03, R spot =39 • , and spot co-longitude spot =178 • (on the side of the primary away from the secondary). The lower panel in Figure 22 plots the radial velocity curves for each model, illustrating that either of the T 2 /T 1 =1 models having q≈0.2 produce acceptable fits to the kinematic data while the nominal best-fitting variable-temperature-ratio model yields a mass ratio that is approximately the inverse of the correct one. At the inclination and mass ratio of the spotted model, the velocity semi-amplitudes require M 1 =1.30 M and M 2 =0.23 M . This best-fitting spotted model produces a line profile (green crosses in Figure 21 ) that matches the primary well in amplitude and radial velocity but slightly underpredicts the velocity amplitude of the secondary at both phases. The peak of the secondary's BF is also slightly smaller than that predicted by the model line profile function. Overall, the spotted model produces a satisfactory match between the model velocity profile and the BF data. By contrast, both of the non-spotted models require substantial third light, l=0.39 and l3=0.44, respectively, which should show up as deficits in the model line profile function compared to the BF data. High-resolution echelle spectra of KIC08913061 Cook & Kobulnicky (2022) show no evidence for third light, effectively ruling out the solutions from the first two models. The lack of a physical explanation for a hot spot on the primary led us to explore the possibility of a cool spot on the secondary. Using two symmetrically placed cool spots with T spot /T 1 =0.69 on the secondary near co-long spot =30 • also produces a satisfactory fit to the light curve, however, the model line profile function still fails to match the secondary peak in the BF in detail and the light curve fit is not as good as when using the single hot spot. We ran MCMC simulations for four free parameters to understand the constraints afforded by the light curve alone. Figure 23 presents the distribution of parameters, showing that the key system parameters are reasonably well constrained and that the correct q is identified. The 16th/50th/84th percentile parameters are cos i=0.027/0.083/0.160, f =0.51/0.84/0.96, log q=−0.824/−0.718/−0.595, and l3=0.14/0.28/0.39. The most probable q=0.19, very similar to the q≈0.17 implied by the velocity data and by the best spotted model. The inclination is high and the mass ratio is extreme. The fillout factor is large, and there is a modest third light contribution. We also ran MCMC simulations for a seven-parameter spotted system with cos i, f , q, l3, T spot /R 1 , R spot , and L spot as free parameters. Figure 24 plots the posterior probabilities. The degeneracy between T spot /T 1 and R spot is evident as a long banana-shaped shaded locus. There is also considerable degeneracy between i and l3 and between q and l3smaller third-light fractions demand smaller mass ratios. Mass ratios between 0.17 and 0.25 are probable. However, the light curve alone is enough to show that extreme mass ratios near q=0.20 are most probable, consistent with the velocities of components measured in the broadening function, albeit with less precision. • . Even with the extra free parameters for the spot, the MCMC simulations are still able to place strong constraints on nearly all of the key parameters. In particular, the spot longitude is well-constrained. KIC08913061 is another long-period, extreme-q system and a probable contact binary. Despite having a light curve symmetric about φ=0.5, the nominal fixed-temperature-ratio contact binary model was not able to produce a good match to the data without recourse to an additional component, modeled here either as a variable temperature ratio T 2 /T 1 =0.86 or using an ad hoc hot spot on the primary. Hence, the real geometrical configuration in this system is ambiguous owing to the need for an extra physical component in the model. The physical origin for a hot spot on the primary opposite the secondary is not obvious, as it would be an unusual location for a mass transfer stream to impact. The lack of agreement between the model line profile function and the BF at the secondary's velocity at both orbital phases indicates a remaining deficiency in the model; a larger fillout factor would improve agreement with the BF but at the cost of larger χ 2 ν in the model light curve. We are unable to entirely reconcile the best-fitting model light curve and the line profile function produced by this model with the BF. KIC08913061 serves as a warning that blind variable-temperature-ratio models fitted to single-band light curves can yield grossly erroneous system parameters. In the case of KIC08913061 (morph=0.76), q and T 2 /T 1 turn out to be degenerate in ways that that lead to an incorrect set of system parameters. The situation becomes more ambiguous when additional physical components are present in the system (e.g., spots). The correct q≈0.2 solution, in fact, does constitute a local minimum in parameter space, but the global minimum with q≈4 and T 2 /T 1 =0.86 is deeper. Without a priori kinematic knowledge of the mass ratio, the correct mass ratio in KIC08913061 is only obtained by constraining the component temperature ratio to T 2 /T 1 ≈1. However, fixing T 2 /T 1 =1 is likely to be overly constraining Figure 24 . Posterior probability distributions for combinations of seven free system parameters for the KIC08913061 system, including a hotspot on the primary. There is considerable degeneracy among i, q, and l3, but the models still place strong constraints on the system parameters. Figure 25 displays the Kepler light curve of the P =1.11 d system KIC09164694, which shows greater complexity than any previous examples. Primary eclipse that is much deeper than secondary eclipse. The maximum near φ=0.25 is significantly brighter than that at φ=0.75. The red curve is constructed using the first nine Fourier components, as labeled in the figure. Some dispersion about the mean light curve is apparent in a small fraction of the time baseline. The minima have flattened bottoms, suggesting an eclipse of both components. Figure 26 presents the eclipse timing residuals, showing a systematic trend that is represented either by a parabolic function implying a period decrease dP/dt=7.7×10 −9 or a periodic function (implying an orbit about a third body) with period P =6.3 yr, amplitude 1.1 minutes, and t 0 =2455676 (time when the contact binary is at superior conjunction). Given the limited time baseline of the Kepler dataset, it is not possible to distinguish between these possibilities, or a third scenario involving quasi-random orbital period modulations arising from magnetic cycles (e.g., Applegate 1992) . We fit the BF with a two-component Gaussian function to measure velocities and relative luminosities, finding it necessary to constrain the velocity width of the fainter component to be 1/3 that of the brighter (140 km s −1 ) 21 , but allowing the centers and normalizations to vary independently. This leads to a radius ratio of R 2 /R 1 0.26, assuming that the temperatures are similar-an assumption that will be relaxed in the analysis that follows. The ratio of velocity semi-amplitudes indicates a mass ratio q 0.20. Because of the BF component blending, the velocities and relative areas are more uncertain than in previous examples where the components are better-separated. The velocity of the fainter component is considerably uncertain at both quadrature phases. Figure 28 shows light curves, residuals, and velocity curves of KIC09164694 for three competing models. The light curve is complex, requiring a great degree of fine-tuning to reproduce in detail. The flattish minima require high inclinations. However, such large inclinations lead to very deep minima-much deeper than observed-a problem that can be remedied by including a modest third-light contribution. The brighter maximum just after φ=0.25 (when the primary component is approaching) compared to just before φ=0.75 requires the introduction of an asymmetry in the system. This can be achieved by placing hot spot on the leading face of the primary or the trailing face of the secondary. We adopt the former, admittedly without physical motivation. Placing the spot on the larger star allows a less extreme T spot /T 1 ratio and a smaller spot size than if the spot were placed on the secondary. Figure 28 shows that the nominal variable-temperature ratio model (magenta dashed curve) yields a reasonable fit (RMS=0.0209) for i=70.5, f =0.84, q=10.1, l3=0.37, and an extreme ratio T 2 /T 1 =0.71. Best-fitting detached models have a much larger RMS and require that both components overfill their Roche lobes, so they are not considered further. Fixed-temperature-ratio models (red dotted curve) produce an RMS that is twice as large, but for very different system parameters: i=89.2, f =0.99, q=0.24, l3=0.36. Only the latter model is consistent with the mass ratio obtained from the kinematic data (q≈0.2). A much better fit (RMS=0.0039) is obtained with the spotted variabletemperature-ratio model (green curve) for i=89.6, f =0.61, q=0.24, l3=0.22, T 2 /T 1 =0.82, T spot /T 1 =1.03, R spot =53 • , and L spot =233 • (roughly on the leading face of the primary). At the cost of three additional free parameters, the reduction in RMS for this model is dramatic. Even so, the best-fitting inclination and mass ratio are nearly identical to the fixed-temperature-ratio model. Reassuringly, the mass ratios of both of these latter models agree with the BF-the spotted model and the T 2 /T 1 =1 model both reproduce the observed velocities of the components at quadrature phases. Including irradiation effects changes the best fits negligibly. We also performed MCMC simulations of the full eight-parameter model (variable temperature ratio plus a spot) to understand the extent to which system parameters can be constrained in the presence of additional physical features in the system. Figure 30 presents the Bayesian distribution. As above, the most probable inclinations and mass ratios are consistent with the kinematic data and the allowed range of other parameters remain consistent with the non-spotted model. There is considerable degeneracy between the spot radius and temperature, but other parameters are well-constrained. Bayesian analysis of KIC09164694 illustrates that the addition of a single physical feature (a spot with three free parameters) greatly improves the model fit and simultaneously provides tight constraints on the principal system parameters and spot parameters. KIC09164694, like many systems in the parent sample, has a complex asymmetric light curve (morph=0.75) that is not well-fit using a standard four-parameter contact binary model with T 2 /T 1 =1. The five-parameter variabletemperature-ratio model yields a best fit with a wrong mass ratio (q=10 versus the correct q≈0.24 measured from the BF), further illustrating the degeneracy between q and T 2 /T 1 revealed above in the case of KIC08913061. Although there is strong evidence for being a contact binary, the large disparity in retrieved component temperatures T 2 /T 1 =0.82 (even with the inclusion of a hot spot on the primary!) seems implausible for a true contact system. We regard this system as an ambiguous geometry based on the present data. The presence of some third light indicated by all models is consistent with the O − C eclipse timing residuals. Although a third component comprising ∼0.2 of the system luminosity is consistent with the BF, but such a stationary component may easily be masked by the brighter primary component. At q≈0.24, KIC09164694 is another extreme-mass-ratio (near-)contact binary on the long-P tail of the period distribution. Figure 30 . Probability distribution for combinations of eight free system parameters for the KIC09164694 system. Even with the additional free parameters entailed by a spot and T2 < T1, the Monte Carlo simulations provide a well-constrained measure of the q that is consistent with the spectroscopic data. Figure 31 displays the Kepler light curve of the P =1.04 d contact system KIC09345838. Like KIC09164694, primary minimum is much deeper than secondary minimum and the φ=0.25 maximum is brighter than the φ=0.75 maximum. The eclipse timing residuals show no systematic variation over the baseline of the Kepler mission. The BF of KIC09345838 is asymmetric with a positive wing at φ=0.32 and a negative wing at φ=0.73, indicating one dominant component and one much less luminous component. The radial velocity of the fainter component is particularly uncertain at both quadrature phases owing to its low amplitude. The component velocities lead to semi-amplitudes of K 1 =32 km s −1 and K 2 =163 km s −1 . These imply a mass ratio near q≈0.20, with considerable uncertainty given the difficulty in measuring the fainter component. Figure 33 displays the light curves and velocity curves of three competing contact models in comparison to the data. The results mirror those seen previously with KIC09164694. Detached models are ruled out, as they all require maximum radii for both components that exceed their respective Roche limits. A variable-temperature-ratio model (magenta dashed curve) provides a reasonable fit (RMS=0.0073) for parameters i=71.5 • , f =0.79, q=1.35 (inconsistent with the kinematic data), l3=0.63, and T 2 /T 1 =0.71. The great disparity between primary and secondary minumum drives this large temperature difference. A fixed-temperature-ratio mode yields a significantly worse fit (RMS=0.0255) for a larger inclination (i=89 • ), f =0.99, a vastly lower mass ratio (q=0.19), and l3=0.50. Only the latter model is consistent with the kinematic data in Figure 32 (q≈0.2. A T 2 /T 1 =1 spotted model provides a superior match to the data (RMS=0.017) with i=75.9 • , f =0.89, q=0.24, l3=0.47, T spot /T 1 =1.06, R spot =60 • , and L spot =190 • (on the primary opposite the secondary). Either of the models with equal-temperature components correctly predicts the mass ratio. The modest l3 in these models is also consistent with the deficit seen in the model BF (green x's) relative to the data, which indicates substantial third light centered near the primary's velocity. Figure 34 shows the posterior parameter probabilities resulting from Monte Carlo retrievals for four free parameters: cos i=0.054/0.149/0.281, f =0.39/0.70/89, log q=−01.08/−0.858/−0.569, and l3=0.14/0.33/0.56. The most probable third-light fractions near l3=0.33 are consistent with the (needed but unmodeled) contribution to the overall line . Three competing contact binary models for KIC09345838 in comparison to the data. As in KIC09164694, the variable-temperature-ratio model yields an incorrect mass ratio, while the fixed-temperature-ratio models yield q≈0.24, in agreement with the kinematic data. profile that would reconcile the theoretical profile with the BF in Figure 32 . Despite some degeneracy between log q and l3, the mass ratio and inclination are well-constrained and the fillout factor is large. . Posterior probability distributions for combinations of four free parameters in the KIC09345838 system. The most probable q of 0.14 provides a good match to the kinematic data. Like KIC09164694 there is some degeneracy between q and l3. Figure 34 shows the posterior parameter probabilities for a seven-parameter T 2 /T 1 =1 spotted model. All parameters except f are well-constrained. The allowed ranges overlap with those from the simpler but poor-fitting four-parameter models. KIC09345838 is another long-period extreme-q system (morph=0.75) where the rather extreme mass ratio is constrained on the basis of the light curve alone, even with the additional model complexity of a stationary spot. However, like KIC08913061 and KIC09164694, the flexibility of the variable-temperature-ratio model leads to an incorrect mass ratio, while a T 2 /T 1 ≈1 model recovers a mass ratio q≈0.2, consistent with the kinematic data. (Smith et al. 2014 ). The parabolic fit under the assumption of a constant period derivative (red curve) fits the data less well than the sinusoidal fit (dashed green curve), which provides an excellent match. The best-fitting sine curve parameters indicate a semi-amplitude of 65.3 min, equivalent to a projected semi-major axis a 1 sin(i) for the contact binary about the system barycenter of 7.8 AU. The best-fitting time of superior conjunction (when the contact binary is farthest away and the third component is nearest) is t 0C =2457400. The best-fitting period is 18.67 years. For an assumed total mass of the inner contact binary of 1.6 M , the minimum implied mass for the unseen third component in a circular orbit would be 3.0 M -sufficiently large and bright for any main-sequence star that it should be the dominant spectral signal in the broadening function (which may be the case, as discussed below). This lower mass limit also rules out white dwarfs and ordinary stellar mass neutron stars as companions. Figure 38 shows the broadening function for two epochs of spectroscopy obtained near phases φ=0.27 and φ=0.74, assuming a quadratic ephemeris with P =0.8784799 d and dP/dt=−3.268×10 −8 and t 0 =2454954.556810. 22 Two components are visible, one much more luminous than the other. The estimated velocity semi-amplitudes are K 1 =30 km s −1 and K 2 =192 km s −1 at a systemic velocity of γ=−43 km s −1 , yielding q≈0.16. The ratio of component areas in BF implies R 2 /R 1 ≈0.4, but this is too large if third light is present. The dominant peak in the BF likely contains a luminous third star, meaning that the velocity amplitude of the primary star in the contact binary is underestimated, making the mass ratio a lower limit. The plotted green x's from the best-fitting light curve model are for q=0.27. The best-fitting variable-temperature-ratio contact model without irradiation (magenta dashed curve) shown in Figure ratio model (RMS=0.0015) is very similar at i=80.7 • , f =0.99, q=0.32, l3=0.63, leading to R 2 /R 1 =0.66, R 1 =2.22 R , R 2 =1.46 R . Both are consistent with the BF results for q and for l3. Models implementing irradiation effects have a similar slightly worse RMS but similar best-fitting parameters. Best detached models have a much larger RMS and are not considered further. Adopting the best model i and q and the observed K 2 =192 km s −1 yields M 1 =0.89 M and M 2 =0.24 M . The model velocity curve in the lower panel of Figure 39 shows the the BF data underestimate the amplitude of the primary, ostensibly because of the bright tertiary that comprises ≈63% of the light and is blended near zero velocity with M 1 . The presence of the third component is unsurprising, given the eclipse timing variations. Figure 38 relative to the model line profile curve. The radial velocity of the primary is also offset from the peak of the BF, in the direction consistent with the third star having a velocity near the system's center of mass. The discrepancy is consistent with the third star providing about 60% of the system light. Nevertheless, on the basis of the light curve alone the principal system parameters are well-constrained, as illustrated by Figure 40 . Compared to the other long-period contact binaries discussed previously, the most probable mass ratio of q=0.31 appears rather modest. KIC09840412 is a clear case of a contact binary (morph=0.79) with nearly equal-temperature components. Variabletemperature-ratio models provide a somewhat better fit and produce a slightly shifted locus of system parameters. Bayesian retrieval of system parameters are consistent with the kinematic data. Figure 41 displays the Kepler light curve of the P =1.38 d system KIC09953894, displaying both both strong ellipsoidal modulation and eclipses. Eclipse timing residuals are <0.2 minutes and show no systematic pattern, indicating that a linear ephemeris is adequate. The light curve morphology parameter of 0.74 is the lowest in our pilot sample and consistent with the possibility of a detached system, indicated also by the V-shaped inflection near primary and secondary minimum. Figure 42 shows the broadening function of KIC09953894 for two epochs of spectroscopy obtained at phases φ=0.23 and φ=0.69. Two well-separated components indicate a primary that is more luminous than the secondary. Component velocities are well-measured and yield K 1 =115 km s −1 , K 2 =141 km s −1 , γ=−32 km s −1 , and q=0.81. The ratio of component radii is near R 2 /R 1 =0.77 if the temperatures are similar. The best model line profile, however, provides a poor match to the data and underpredicts the mass ratio and the velocity semi-amplitude. Figure 43 presents the model light curves and velocity curves for three competing models in comparison to the KIC09953894 data. The nominal fixed-temperature-ratio model (magenta curve) provides a reasonable fit (RMS=0.0028) for i=86.7 • , f =0.77, q=0.96, and l3=0.79, but the residuals are significant and systematic particularly surrounding φ=0.5, as the contact model is unable to reproduce the inflection indicative of eclipses. The 16th/50th/84th percentile ranges are especially broad (q=0.78-1.17). The model does, however, approximately reproduce the velocity semi-amplitudes of the components in the lower panel. A detached model (green curves) provides a much better fit (RMS=0.0014) for a much lower inclination i=60 • , similar q=0.90, much lower l3=0.00, and T 2 /T 1 =1.06 (a hotter secondary). The middle panel shows that the residuals for this model are dominated by imperfect agreement near φ=0.5. Including irradiation and reflected light (blue curve) reduces the residuals near φ=0.5 (RMS=0.0011) but does not eliminate them, suggesting either additional physical effects in the system or an imperfect implementation of irradiation in the PHOEBE models. The resulting parameters are similar to the detached model except that q=0.97, closer to unity than allowed by the kinematic data (q=0.81). However, q=0.81 is within the ±1σ uncertainties in all three models. In both of the detached models the best-fitting parameters yield R 1 /R 1max >0.99, indicating a primary close to overflowing. Both the detached models underpredict the velocity semi-amplitudes; this is a consequence of fixing M 1 =1 M in the models, while the kinematic data dictate M 1 =1.3 M for i=60 • . With an indicated primary temperature of 7295 K, KIC09953894 appears to be a relatively hot system probably containing two late-F stars, which would be consistent with the masses derived from the spectra, assuming a main-sequence temperature-mass relation. The blue curve shows the best-fitting detached model with irradiation. This model provides a slightly better fit than the detached model without irradiation effects, but the residuals near phases φ=[0.0, 0.5] are still systematic. KIC09953894 serves as a warning that for detached systems light curves alone are insufficient to obtain an unique solution; kinematic data are required to fix either q or M 1 . Even so, the Bayesian confidence interval limits q to values that include the spectroscopically measured mass ratio. Figure 44 shows the MCMC distribution of system parameters for the KIC09953894 system for the detached model. Percentile ranges are cos i=0.482/0.500/0.517, log q=−0.105/−0.054/−0.018, l3=0.00/0.02/0.04, and T 2 /T 1 =1.02/1.03/1.04. By comparison the most probable parameters for a worse-fitting contact model are cos i=0.065/0.125/0.205 (much larger than the detached model), f =0.46/0.67/0.85 (not applicable in the detached model) log q=−0.106/−0.008/0.067 (similar to the detached model), and l3=0.72/0.76/0.79 (much larger than the detached model). While the mass ratio distribution encompasses the correct q∼0.81, the ranges for the preferred detached model are very different than the contact model despite a similar RMS. The large third light fraction retrieved from the contact model is inconsistent with the l3≈0 indicated by the detached models and by the BF. The secondary is slightly hotter than the primary (T 2 /T 1 =1.03) and the primary is close to overflowing (R 1 /R 1max =0.99) but the secondary is not (R 2 /R 2max =0.8), yielding large ellipsoidal modulations of the light curve and shallow eclipses at both quadrature phases. 2 .4 8 2 .5 2 2 .5 6 2 .6 0 R1 0 .6 0 .7 0 .8 0 .9 R2/R1 0 .9 8 7 0 .9 9 0 0 .9 9 3 0 .9 9 6 0 .9 9 9 R1/R1max 0 .6 0 .7 0 .8 0 .9 R2/R2max Figure 44 . Bayesian probability distribution of KIC09953894 system parameters when modeled as a detached configuration and a primary mass of M1=1 M is assumed. The data and models indicate that KIC09953894 is a detached system with two inflated components and a slightly hotter but smaller secondary. Contact models can approximate the light curve and yield a mass ratio broadly consistent with detached models and spectral data, but the contact models erroneously predict large third light and larger inclinations. Including irradiation and reflected light in the models improves the agreement with the light curve for this T 2 /T 1 =1.03 system. Without a spectrosopcially determined primary mass or mass ratio, both detached models (with and without irradiation effects) underpredict the velocity semi-amplitudes but do an excellent job of reproducing the light curve. Figure 45 displays the Kepler light curve of the P =0.56 d contact system KIC10292413, the shortest period among our pilot sample. The primary minimum is very slightly deeper than secondary minimum. This target also has a close visual companion within 2 which may affect the measured flux of the Kepler data since it lies within one detector pixel. Dilution of a light curve by a constant third light makes the maxima fainter and the minima brighter, reducing the overall amplitude of modulation, in agreement with the fairly small amplitude in Figure 45 . The Pan-STARRS (Kaiser et al. 2010) images of this target in Figure 46 show it to be a visual double consisting of approximately equal brightness components at a position angle of 120 • at a 2. 3 separation. We spectroscopically determined that the indicated star to the southeast is the contact binary, identified as Gaia DR2 ID 2086341391231975168 (G=14.6 mag). We aligned the spectrograph slit angle to cover only this star and minimize contributions from the neighbor to the northwest (Gaia DR2 ID 2086341391231974656; G=15.1 mag). Figure 47 presents the eclipse timing residuals of KIC10292413, which show some evidence of a decreasing period with dP/dt=−1.06×10 −9 (red parabola). This is much smaller than in previous examples and is at the limit of what can be detected in the data. A sine function with a period of 10 yr is also a good fit to these data. The amplitude of 0.78 minutes corresponds to a minimum projected light crossing distance of 0.1 AU for the orbit of the contact binary about a common center of mass with a hypothetical third body. A very low-i orbit about a M 3 ≈0.5 M star could plausibly explain these O − C variations. We conclude that the evidence for eclipse timing variations is strong in this object, but we are unable to distinguish between a constant period derivative and a periodic O − C variations that would indicate a third body in the system. In any case, the visual companion at 2. 3 separation cannot constitute the third body responsible for the eclipse timing variations, as it has a very different Gaia parallax and must be physically unrelated. Figure 48 shows the broadening function of KIC10292413 for two epochs of spectroscopy obtained at phases φ=0.23 and φ=0.73 (adopting the best linear ephemeris). Judicious slit placement allowed us to exclude light from the close visual companion. Owing to the faintness of this source (G=14.6 mag) the signal-to-noise ratio is lower than for other targets. Most of the information in the broadening function comes from the Hα line which is the strongest in the spectrum. There is no reported T eff in the Kepler Input Catalog, so we adopt 6200 K, a value consistent with the strong Hα line. There are two peaks in the (rather noisy) BF, allowing us to estimate component velocities and a mass ratio near q≈0.37, with considerable uncertainty. The estimated velocity semi-amplitudes are K 1 =70 km s −1 and K 2 =189 km s −1 at a systemic velocity of γ=−28 km s −1 . The ratio of component areas in BF implies R 2 /R 1 ≈0.71. However, the presence of third light likely skews the velocities of one or both components, rendering the kinematic measurements significantly uncertain. Figure 49 presents the best-fitting fixed-temperature-ratio contact binary light curve and velocity curve for KIC10292413. A excellent fit (RMS=0.00010) with no systematic residuals is achieved with a contact binary system using i=53.7 • , f =0.88, q=0.53 (somewhat larger than the value inferred from the BF), and l3=0.59, leading to R 2 /R 1 =0.79, R 1 =1.68 R , and R 2 =1.32 R . At this inclination the velocity semi-amplitudes imply component masses M 1 =1.34 M and M 2 =0.71 M . Detached configurations produce a much larger RMS and are not considered further. Models employing a variable temperature ratio provide an essentially identical solution, supporting equal-temperature components. The model line profile (green x's in Figure 48 ) matches the BF reasonably, despite the lower SNR of the spectra. Yet, the deficit between the BF and the model line profile suggests a third component at a different radial velocity in the two spectra obtained more than a year apart. Two sources of third light are possible. The visual 2 companion to the northwest constitutes 38% of the combined system light, as it falls within one Kepler pixel. The O − C variations also suggest third light within the KIC10292413 system. Together, these could easily account for the l3≈0.59 suggested by the best model. Despite the large light curve morphology parameter of 0.95 that might suggest a detached ellipsoidal variable, KIC10292413 appears to be a genuine contact system with a substantial third light contribution. KIC10292413 may be a tight triple system with a velocity-variable third component helping to generate the irregular shape of the BF in Figure 48 . Figure 51 displays the folded Kepler light curve of the P =0.999 d contact system KIC11097678. The secondary minimum is not as deep as the primary minimum, and the minima are flat-bottomed, indicating an eclipsing system seen at high inclination angle. Figure 52 presents the eclipse timing residuals of KIC11097678, which show some evidence of an increasing period with dP/dt=2.4×10 −9 . Alternatively, the O − C data can be fit with a sine function having amplitude A=0.96 min, P =8 yr, T 0 =2476356. The short duration of the Kepler dataset precludes more robust conclusions regarding the eclipse timing variations. Unfortunately, the ZT F data are not helpful because of the limited phase range sampled on this P ≈1.00 d system. Figure 53 shows the broadening function of KIC11097678 for two epochs of spectroscopy obtained at phases φ=0.23 and φ=0.79. There are two clear peaks in the BF. The velocity of the fainter component at φ=0.23 is not wellconstrained on account of the smaller signal of the secondary in the BF. The ratio of areas of the two Gaussian components implies a radius ratio R 2 /R 1 =0.35. The radial velocities of the contact components give velocity semiamplitudes K 1 =8 km s −1 and K 2 =252 km s −1 with γ=-74 km s −1 , implying a mass ratio near q∼0.03. The green x's depict the model line profile function of the q=0.10 model that provides the best fit to the light curve, but we find that a model with an even lower q would improve the fit to the BF. The velocity of the secondary is not well measured at the first quadrature phase, but an extreme mass ratio is clearly indicated by the velocities. At this phase the BF shows a substantial deficit relative to the model, consistent with a real physical effect that decreases the flux from the secondary. If there were a third component blended with the primary near the systemic velocity, our measurement of K 1 and q would be lower limits. The Pan-STARRS image data show a visual companion 1. 5 to the west, 3.3 mag fainter in the Gaia photometric system, implying a third-light contribution of about 5% within the Kepler pixel. Such a faint companion does not contribute measurably to the spectra or show up in the BF. It has a substantially larger parallax and cannot be physically associated with KIC11097678. Figure 54 shows a light curve and velocity curves for KIC11097678 along with a best-fitting (RMS=0.0024) variabletemperature-ratio contact model (magenta curve) having i=84.5 • , f =0.92, q=0.10 (somewhat larger than the q≈0.03 measured from the BF), l3=0.31, and T 2 /T 1 =1.02 (a slightly hotter secondary improve the fit. However, the fixed-temperature-ratio contact model provides a similarly good fit (RMS=0.0024) for i=65.9 • , f =0.62, q=0.06, and l3=0.01. Competing detached models imply R 1 /R 1max and R 2 /R 2max 0.99 and are not viable, having RMS at least three times larger. The best model component radii are R 1 =3.36 R and R 2 =1.05 R , yielding R 2 /R 1 =0.31, very similar to that inferred above from the BF. Component masses approximately M 1 ∼1.93 and M 2 ∼0.19 using the velocity data and best model i. This PHOEBE model produces a line profile function (green crosses in the BF plot of Figure 53 ) that is in excellent agreement with the BF. (in fact, the significant third light implied by the variable-temperature-ratio model (with or without irradiation) cannot be accommodated in the BF, lending credence to the T 2 /T 1 =1 model!) The agreement with the BF is good in the lower panel at φ=0.79 except that the secondary peak has a larger amplitude than the data. In the upper panel at φ=0.23, the primary's BF peak is shifted to positive velocities relative to the model and the secondary appears too faint in the BF compared to the model profile, echoing the deficit noted previously for KIC04999357 and KIC10292413. This difference exists regardless of the range of wavelengths used to construct the broadening function. We conclude that it is not a consequence of Hα emission or some other portion of the spectrum that varies with orbital phase. 23 We find that the model line profile at φ=0.23 can better reproduce the BF if we add a large cool spot (T spot /T 2 =0.8; or indeed any phenomenon that suppresses the signal from the secondary) on the trailing face of the secondary (co-longitude ∼90 • ). Such a spot configuration would not alter the good agreement between the model and BF at φ=0.79. Large transient spots are known to exist on contact binaries (e.g., Binnendijk 1970; Maceroni et al. 1994; Kouzuma 2019) . Although such large spots would be visible in the light cure, it is possible they were not present (or present only for a brief duration) during the 2009-2013 Kepler epochs. Figure 55 shows the relative probabilities resulting from Monte Carlo simulations using the four-parameter fixed-temperature-ratio models. The percentile parameters are cos i=0.195/0.225/0.251, f =0.47/0.58/0.67, log q=−1.19/−1.16/−1.12, and l3=0.01/0.05/0.10, implying strong constraints on system parameters. These parameters are consistent with those of the best fits shown in Figure 54 , indicating very small q and small third-light fraction. The most probable l3 is consistent with the estimated 5% light contribution from the identified 1. 5 visual companion (GDR2 2086752574218680704). Contours in Figure 55 are helpful in illustrating the degeneracy between log q and l3, here in the sense that larger third light would lead to less extreme mass ratios. Including irradiation effects in the models changes the RMS negligibly but substantially shifts the locus of most probable parameters toward larger inclinations, larger fillout factors, and larger third-light fractions. Figure 56 shows the posterior distribution of parameters from MCMC simulations where irradiation effects are included. The most probable parameters are significantly different than in Figure 55 where irradiation is disabled: cos i=0.039/0.105/0.183, f =0.69/0.83/0.93, log q=−1.12/−1.05/−0.99, and l3=0.16/0.24/0.32. The non-negligible l3 found when irradiation is included is inconsistent with the BF data in Figure 53 and with our BF analysis of high-resolution echelle spectroscopy (to be presented in Cook & Kobulnicky 2022) which limits the third-light fraction to <0.01. The O − C variations may, then, result from a very faint low-mass companion or from magnetic activity. The substantial third-light fractions implied by irradiation-enabled models cast doubt upon the blind application of this particular implementation, at least in this particular system. Our analysis of KIC11097678 suggests that caution is appropriate when invoking irradiation physics, but we are not able to speculate as to the reason for this inconsistency. Nevertheless, the most probable mass ratios are still extreme for any model-based Bayesian analysis-near nominal limits for onset of the Darwin instability. Despite the large morph=0.90, KIC11097678 is a clear case of a contact system. In the pilot sample, spectroscopic data provided radial velocities of (near-)contact binary stars that permitted the measurement of component masses when combined with orbital inclinations inferred from the light curves. The systems' mass ratios were already well-constrained by the light curves alone, even in the presence of additional features modeled as spots on the stars (KIC08913061, KIC09164694, KIC09345838). The spectroscopic data, nevertheless, provided confirmation of the mass ratios and the third-light fractions retrieved from the light curve analyses. Sometimes the third component could be seen directly as a distinct peak in the BF. Most often the third component was blended in velocity space with the brighter and more massive primary component, adversely affecting the measurement of the primary star's velocity and systematically decreasing the mass ratio inferred from the BF. Notably, all but one of the ten pilot study systems show some evidence for a third component either from the Monte Carlo light curve analysis, the O − C variations, the BF, or some combination of the three. In all but one of the systems, the evidence for a third-light contribution is compelling. Light curve shapes are insensitive to the mass of the primary, so kinematic data are required if component masses are to be measured. We performed PHOEBE+MCMC simulations of the KIC09345838 system (which required a hot spot on the primary star to achieve good fits to the mean light curve) allowing the mass of the primary, M 1 to vary in addition to the seven parameters modeled in Figure 34 . The posterior probability distributions of all parameters are nearly identical to when M 1 is fixed, and the uncertainties are as well. This is expected and merely serves to confirm that the light curve morphology is sensitive to the mass ratio and not to the overall system mass. Table 2 provides a short synopsis of the best-fitting model parameters and spectroscopically determined parameters for the (near-)contact systems in our pilot study. Column 1 is the KIC ID, columns 2-6 give the i, f , q, l3, and T 2 /T 1 from the best-fitting contact binary model, columns 7-10 give the BJD, observed orbital phase, and component Heliocentric radial velocities for the first spectroscopic observation near φ=0.25, columns 11-14 give the BJD, orbital phase, and component Heliocentric radial velocities for the second spectroscopic observation near φ=0.75, column 15 lists the computed systemic velocity, and column 16 gives the spectroscopic mass ratio, q spec , if measurable from the BF. Table footnotes provide additional insights on each system gleaned from the spectroscopic analyses. Five systems from the pilot study (KIC04999357, KIC06844489, KIC09840412, KIC10292413, and KIC11097678) have contact models that provide an excellent match to the light curves and have RMS considerably less than the best detached models. Variable-temperature-ratio models provide slightly better fits and have 0.980.70 indicative of contact configuration or ellipsoidal variations from at least one of the aspherical components. For 45 systems we assigned a different reference time T 0 shifted by half an orbital period compared to the one tabulated in Kirk et al. (2016) in order to place the deeper eclipse at phase zero. We removed from the sample KIC02831097, as it is a pulsating RR Lyrae type star described by Sódor et al. (2017) . We also removed 23 systems owing to their extremely low level of photometric variability where the amplitude of modulation is less than about 10 times the RMS in each phase bin. Such systems are likely to be low-i and/or high-l3 systems, and the data do not yield useful constraints on system parameters. We further removed KIC07733540 which showed large O − C variations and abnormally large photometric variability. This left 783 systems to model. We performed MCMC simulations of all 783 systems, modeling each as both a T 2 ≈T 1 (constrained temperature ratio) contact binary using five free parameters and as a detached system using six free parameters. Both sets of models include irradiation effects. Although contact binaries can exhibit small temperature differences between the components, the degeneracies documented with the pilot sample using the variable-temperature-ratio contact model motivated us to abandon this unconstrained temperature ratio option as dangerously unreliable. Photometric data from a second bandpass could, presumably, break this degeneracy and retrieve all of the system parameters with better fidelity. No single criterion can reliably discriminate between contact and detached geometries. We employ three model-based metrics that provide approximate classification as probable contact binaries (C), probable detached systems (D), or ambiguous cases (A). Our classification scheme (based on lessons learned from the pilot sample) utilizes the RMS from the best fitting contact versus detached models (i.e., RMS con versus RMS det ), the temperature ratio T 2 /T 1 from the best-fitting contact model, and the fraction of the components' Roche lobes filled from the best-fitting detached model (i.e., R 1 /R 1max and R 2 /R 2max ). • (C) Contact-178 systems: RMS con 0.95). This last criterion requires at least one of the components to fill its Roche lobe to at least 95% of the maximum value in the best-fitting detached model. RMS det 0.7) already select against most eclipsing detached systems but do allow ellipsoidal variables which tend to have morph>0.9 (Kirk et al. 2016) . Probable contact systems congregate preferentially at short periods. Ambiguous systems are found at all morphology parameters but cluster strongly at the shortest periods. Long-period systems P >0.8 d are preferentially detached systems which populate the upper right corner at morph>0.9 as expected for ellipsoidal variables. Very few contact systems lie at morph<0.75. Ambiguous systems concentrate at short periods and span the full range of morph but show a concentration near morph∼0.75. Many ambiguous systems receive their designation for having poor fits to either model, generally because the light curve is asymmetric. We suspect that asymmetric light curves end up with morph≈0.75 designations as a consequence of the information loss in downprojecting a higher dimensional manifold to a single dimension (see Kirk et al. 2016, and descriptions therein Figure 57 . Light curve morphology versus orbital period for the probable contact binaries, detached binaries, and ambiguous cases, as denoted by symbol type. Figure 58 presents histograms of the uncertainties on each fitted parameter, characterized as 0.5 times the difference between the 84th and 16th percentile value, i.e., 0.5∆cos(i) 84−16 , for systems modeled as contact binaries. These values constitute approximate 1σ uncertainties as long as the distributions are nearly Gaussian (they are often not). The gray shaded histogram includes all 669 C+A systems. Systems classified as contact (blue), detached (red), or ambiguous (black) are plotted separately by color. The median 0.5∆cos(i) 84−16 is 0.045 (upper-left panel) , and all values are <0.22, indicating that the MCMC simulations generally provide strong constraints on inclination. This is especially true for contact (C) systems which have smaller median uncertainties than the other subsets. The median 0.5∆f 84−16 is 0.15 (upper middle panel) but smaller for the contact subset. The median 0.5∆log(q) 84−16 is 0.11 (upper right panel), with a few cases as large as 0.38, indicating the the probabilistic constraints on log q are strong in most instances-again, much smaller for the contact subset. Detached systems have larger uncertainties, as expected, given that the detached geometry light curve is not sensitive to q. The median third light uncertainties (lower left panel) are 0.5∆l3 84−16 =0.08. The median T 2 /T 1 (lower middle panel) is 0.009, indicating that the temperature ratio is constrained to better than 1% in most cases, although the artificial hard bounds imposed on T 2 /T 1 at [0.95, 1.05] mask the true extent of possible variations. The lower right panel presents the distribution of the RMS between the mean of the data and the best-fitting contact binary model, where the median RMS is 0.0022. This is seven times larger than the nominal Kepler photometric uncertainty of 0.0003, affirming that the data contain the signatures of real physical features not present in the nominal PHOEBE models. However, the distribution peaks sharply at small RMS for all subsets, meaning that a large fraction of the sample exhibit reasonable agreement with the models. The shaded region at right marks the upper limit of RMS=0.005 defining one criterion for the ambiguous designation. Figure 59 plots the Bayesian uncertainties versus cos i, with contact/detached/ambiguous systems represented as blue squares/red triangles/grey dots, respectively. The top panel show that the typical uncertainty on cos i is smallest among contact systems and that the uncertainty decreases toward large cos i (small i; face-on systems). This may be understood as a consequence of all the projected quantities in a binary system (e.g., linear size, velocity) scaling as sin(i), which changes most rapidly at i≈0 • , allowing for better discrimination between models with similar inclination. The uncertainties on f are similar between C and A systems, and there is a modest trend with cos i in the sense that f is better constrained at large inclinations (small cos i). This is where the eclipsing geometry creates light curves laden with information about the system parameters. Constraints on the mass ratio (third panel) are much stronger for C than for A or D systems, unsurprisingly, regardless of inclination. The mass ratio in contact systems becomes less well constrained at large cos i. A similar pattern is evident in the lower panel where the typical uncertainty on third-light fraction grows with cos i. Constraints are similar for all subsets. Taken together, these panels indicate that the best-constrained systems will be those at large i where the (near-)eclipsing geometry provides the mostly highly structured light curves containing information on the principal system parameters. Figure 60 displays histograms of best-fitting model parameters for all 669 systems modeled as contact binaries or ambiguous cases. The upper left panel shows that the distribution of cos i is strongly peaked at small and large values when all C+A systems (gray shaded histogram) are considered. The majority of these are ambiguous systems which have ill-fitting models (RMS>0.005) or large temperature differences between the stars (T 2 /T 1 <0.96 or >1.04). The distribution becomes nearly flat, as expected for random viewing angles, when only probable contact systems are considered (dashed blue histogram) but there is still a statistically significant excess at inclinations cos i≈0.3 (i≈70 • ). Investigation shows that some of these systems have light curve asymmetries and RMS values near the 0.005 cutoff. A more stringent RMS<0.004 cutoff significantly reduces this excess, which we tentatively attribute to the effects of starspots. Furthermore, spots at longitudes 0 • or 180 • in moderate-inclination (cos i≈0.6) systems will change the ratio of primary:secondary minima, displacing them toward more edge-on orientations (smaller), further contributing to the apparent excess of cos i≈0.3 systems. Spots appear to be abundant in contact systems, but a full exploration of spot properties is beyond the scope of this work. It is also noteworthy in Figure 60 that the distribution of probable detached systems-which have poorly measured inclinations, in general-peaks in the cos i≈0.5 regime. At these viewing angles detached but nearly overflowing stars undergo grazing eclipses are most likely to be identified as they produce distinctive V-shaped inflections in the light curve. Detached systems with very low inclinations will exhibit ellipsoidal modulations and appear much like contact systems and may be incorrectly identified as such. As a case in point, KIC09953894 from the pilot study has an indicated i≈87 • when modeled as a contact system but i≈59 • when correctly modeled as a detached system. At the other extreme, detached systems with very large (eclipsing) inclinations have already been removed from the sample by virtue of the initial selection criterion morph>0.7, thereby producing the apparent deficit of cos i≈0 detached systems. The left middle panel of Figure 60 shows that the distribution of fillout factor is slowly declining with f contact systems. There is a small excess near f =1 for contact systems that would indicate a preference for Roche lobes filled near the maximum volume. The upper right panel (log q 50 ) shows that mass ratio distribution for both the C+A and C subsets is similar, evincing two peaks on either side of zero. As discussed in connection with Figure 2 , the log q<0 peak contains systems at inclinations i 50 • while the systems having log q>0 are those with i 50 • . The lower left panel shows that contact binaries have appreciable third-light contributions. The distribution is flat or slowly falling out to l3=0.8 where few systems are found. All subsets show a preference for considerable third light. The distribution demonstrates that 77% of Kepler probable contact binaries have l3>0.15, indicating that luminous triple systems are common in close binaries-much more common than among the general population where this fraction is about 40% among short-period binaries but much lower for longer-period systems (Tokovinin 1997) . Our finding for contact binaries is generally consistent with the 59±8% tertiary fraction discovered among a spectroscopically studied sample of close binaries (Pribulla & Rucinski 2006) . Our result may be slightly biased by our initial removal of the 23 systems having very low levels of photometric variability, some of which are expected to have large l3, but some of these 23 will also be low-i (high-cos i) systems as well. The lower middle panel of Figure 60 shows that a large fraction of the C+A subset is best fit by models having significant temperature differences, T 2 /T 1 <0.96. This is not surprising as this subset contains unidentified detached systems plus a large portion of ill-fitting light curves with strong asymmetries. The probable contact systems, by our adopted criteria, have 0.96≤T 2 /T 1 ≤1.04, and the distribution is nearly flat across the allowed range, with a preference for T 2 /T 1 <1. The distribution of absolute mass ratio in the lower right panel peaks near q a ≈0.4, indicating that the majority of close binaries have mass ratios far from unity, both for the C+A sample and the more restrictive C subset. There are very few contact systems having q a >0.8. Figure 61 plots each of the principle contact binary parameters against one another, with grey points representing the A subset and blue squares marking the probable contact subset. In cases where the mass ratios has been inverted to obtain q a we have also inverted T 2 /T 1 . The pileup of points at T 2 /T 1 =0.95 or 1.05 along the bottom/top in the lower row of panels illustrates the abundance of either detached systems or those having ill-fitting light curves. The distribution of contact binaries shows no strong correlations or patterns across each parameter pair, with a few exceptions. In the q a versus cos i panel there is a paucity of contact systems having moderate mass ratios −0.460 edge-on systems). A substantial population of ambiguous systems do populate this regime. At these large inclinations the constraints on log q are quite strong (Figure 59 ), so we speculate that spots or other physical phenomena not present in the models force the best-fitting solutions toward larger cos i and smaller log(q a ) or cause true contacts to receive an ambiguous designation. The distribution of log(q a ) versus l3 shows that extreme-q a systems are found across the full range of l3, but less commonly when third light is large. This further underscores the warning that system parameters become less certain in the presence of substantial third light. The dependence of mass ratio on orbital period provides the most direct test of evolutionary theories of contact binary systems. The distribution of mass ratio in short-period P <100 d detached systems with ≈1 M primaries (progenitors of contact systems) is observed to be nearly uniform (Niu et al. 2020; Moe & Di Stefano 2017) , so deviations from a uniform distribution in contact systems must be a consequence of evolution. The most commonly cited modification that is expected is a lower limit on mass ratio due to the Darwin instability (Darwin 1893) , generally taken to be a fixed value around 0.09 (Rasio 1995) or ≈0.07 (Arbutina 2009; Li & Zhang 2006) . The evolutionary models of Molnar et al. (2019) ; Molnar (2022) yield orbital period dependence of this minimum, along with period dependences of a maximum q a and of typical q a values. The present data set is well suited to this analysis having 1) the largest unbiased sample analyzed in a uniform fashion, 2) no period dependence in the selection of systems, and 3) no period dependence in any systematic biases in mass ratio determination. In this section we first present the data and then discuss them in terms of the Molnar (2022) evolutionary models, which will be more completely described in that later work. Figure 62 plots the absolute mass ratio q a versus period for the probable contact systems and ambiguous systems modeled as contact binaries. We focus our analysis on probable contact systems represented by the blue squares. The majority of the ambiguous systems, represented by small black points, are probable contact systems as well and can be considered for confirmation of trends, although with more spread in values. Symbols with vertical error bars depict the 50th-percentile and 84th-16th percentile ranges. Each of the three panels shows a different window of mass ratios and periods, allowing the full dynamic range to be conveyed. The upper left panel displays the full range in both parameters, confirming the well-known concentration of contact binary periods in the range 0.25-0.5 days. The mass ratio in this period range is revealed to concentrate in the range 0.2-0.8. There are very few systems with q a > 0.8 at any period. At periods longer than P ≈ 0.5 d the median mass ratio of blue squares shifts toward more extreme values. At P > 0.5 d nearly all of the systems have q a < 0.3. The lower left panel better shows the distribution of mass ratios at longer periods. No systems are found at q a < 0.08. This lower bound becomes more clear when some additional extreme systems from the literature are added. AW UMa, for instance, is a well-studied short-period (P =0.43 d) system with a spectroscopically determined mass ratio of q =0.099 (Rucinski 2015) . A magenta star marks its position in Figure 62 , lying near the lower limit in the q a versus P plane. An additional star marks KR Com (Rucinski et al. 2002) , another spectroscopically confirmed extreme-q system that lies on this boundary. While a few other extreme q a =0.05-0.10 candidates exist in the literature, they are based on photometric measurements only, which carry larger uncertainties on account of the (often unknown) third-light contributions. The lower left panel also shows that long-period contact systems are rare and have q a >0.1 when P ≥1 d, avoiding the most extreme mass ratio regime. The upper right panel shows the short-period portion of the range at higher resolution. It shows the few systems with q a > 0.8 are among the shortest period systems (P 0.3 d). Figure 62 . Colored curves depict evolutionary tracks for contact systems with initial primary mass M 1 =1 M and initial mass ratios q ini =0.1-0.9, color coded and labeled (based on the mass and angular momentum conserving models of Molnar et al. 2019; Molnar 2022) . The stellar contact components begin as evolved single stars at an age of 1 Gyr, based on the >0.6-3 Gyr required for Lidov-Kozai cycles and magnetic braking to bring close detached systems into contact Stepien (1995) ; Chen et al. (2016) ; Hwang & Zakamska (2020) . Increasing or decreasing the primary mass shifts the model tracks to the right or left, respectively, by a small amount. Diamonds of increasing size along each track mark the progression of 2 Gyr intervals from 1 Gyr to 11 Gyr, from upper left toward lower right, as the system transfers mass to the more massive component. The dashed curve defines a lower limit on q and an upper limit on P at which the suite of all model systems (including those of larger and smaller M 1 not depicted) reach the critical threshold for the onset of the Darwin instability when the rotational angular momentum of the components exceeds one-third the orbital angular momentum. The short section of dotted curve in the lower left designates a limit for a small subset of initially very extreme mass ratio systems where the conditions for merger are satisfied upon initial contact. The distribution of systems in Figure 63 , both contact and ambiguous, stands in good agreement with the model predictions. Short-period systems with q≈0.8-1 at the upper left are rare, consistent with the expectation of rapid mass transfer during the first 1 Gyr in such systems, the result of a previously unnamed instability discussed further in Molnar (2022) . The short-period bound along the left edge seen in the data is reflected in the models. The small number of points that lie to the left of model tracks can be reproduced by systems having a slightly smaller total mass or by a slightly younger system age. The upper right portion of the plot is devoid of points, except for a few ambiguous systems that are likely to be detached systems. A small concentration of contact systems near 0.4 d and q a =0.6 lying to the right of the M tot =1.9 M model track may be more slightly more massive systems. No contact systems lie below the dashed line demarcating the locus of stellar mergers; KIC03442006 at P =1.47 d comes the closest, but its 1σ uncertainties straddle the limit. One ambiguous system, KIC09637265, lies below this limit at P =1.85 d, a location where detached systems are permitted. Overall, the bulk of the data points lie at short periods and modest mass ratios-locations where the models predict contact binaries spend the majority of their lifetime. Figure 63 implies several predictions that can, in principle, be tested by detailed studies of individual systems or groups of systems. Systems near the top of the Figure at q≈1 are expected to be young and rapidly evolving in q, perhaps showing signatures of vigorous mass transfer. Long-period contact systems are expected to exhibit extreme mass ratios. Any long-period systems found below the predicted limit are expected to be detached systems. Shortperiod contact binaries near the left edge of the model tracks are expected to be young, achieving contact only recently. The lower left corner of the q versus P should be populated by systems with small total mass and the upper right populated by large-M tot systems. A better knowledge of outliers and systems at the periphery of the model tracks will help inform the boundary conditions for additional modeling efforts and perhaps identify different contact binary formation/evolution channels than the one envisioned in this set of models. Figure 63 . Mass ratio versus orbital period as in Figure 62 . The colored tracks show the evolution of contact binaries-from upper left toward lower right-of contact systems with initial mass ratios qini=0.1-0.9, as labeled. Diamonds mark equal time intervals every 2 Gyr. Decreasing/Increasing the total system mass shifts models slightly to the left/right, as indicated by arrows. Figure 64 compares 50th-percentile parameters as inferred from the detached model with the 50th-percentile parameters from the T 2 ≈T 1 contact model versus log P for the 114 probable detached systems. Unsurprisingly, there are more detached systems at long periods than short periods. The red line in each panel shows the linear fit to the data. The upper panels shows the difference between most probable cos i values, where the RMS dispersion is 0.25. At short periods the differences are symmetric about zero, but at long periods they becomes systematically negative; contact models yield systematically larger cos i (smaller inclinations) than detached models in long-period systems that have detached configurations. This means that detached systems mistakenly modeled as contact systems will have cos i values that are too large (i too small). The second panel plots the difference in log q between detached and contact models. The dispersion is large (0.78) but not systematic, primarily confirming an earlier conclusion that mass ratio is not well constrained for detached systems. The third panel shows the difference in third light fractions, which is consistently negative at all periods with a mean of −0.22. Contact models yield systematically larger third light fractions than detached models; detached systems mistakenly modeled as contact systems will yield l3 values that are too large. The fourth panel plots the ratio of the RMS of the best-fitting detached/contact models; this ratios ranges from near unity to 0.1 across the period range, indicating that sometimes the models are comparable and sometimes the detached model provides a much-superior fit to the light curve. The RMS alone is not a reliable metric to distinguish between contact and detached geometries. The bottom panel shows the average fraction of the components' Roche lobes that are filled in the best-fitting detached model. This derived parameter is akin to a fillout factor for detached systems. At the shortest periods this fraction is near unity, indicating that both stars nearly overflow. At longer periods this ratio increasingly departs from unity, indicating that one or both stars are not close to overflowing. Such long-period systems produce the classical ellipsoidal light curve modulations. Figure 57 has already shown that most of the long-period systems are detached and have morph>0.9. Figure 64 illustrates that the best detached solutions can be very different from the best contact solutions, highlighting the necessity of making a correct geometrical identification that dictates the light curve modeling strategy. Figure 64 . Differences or ratios of 50th-percentile parameters comparing detached models and the T2≈T1 contact models for the 114 probable detached binaries. Parameters from best detached model can be very different from the best contact model. We have presented an analysis of 783 P <2 d close binaries using MCMC methods in conjunction with the state-ofthe-art PHOEBE software to model Kepler light curves. The reliability of the Bayesian retrievals of system parameters is predicated on the success of the models demonstrated using a pilot sample of ten systems for which phase-resolved spectroscopy provides independent measurements of the systems' mass ratios and estimates of third-light fractions. The level of intrinsic variability in all of these systems far exceeds the exquisite photometric precision of the Kepler data, indicating the ubiquity of intrinsic variability from phenomena that may include pulsations, spots, and active regions that vary both in time and location upon the stellar surfaces. Accordingly, the best-fitting and 16th/50th/84thpercentile parameters presented reflect the mean properties as averaged over the three-year baseline of the Kepler dataset. Our detailed photometric and spectroscopic analysis of the ten-object pilot sample provides a representative tour of the types of systems and physical phenomena comprising the Kepler close binaries: contact and detached systems with a variety of fillout factors, spots, temperature differences, O − C variations, and third components. 1. PHOEBE models plus Markov-Chain Monte Carlo analyses of single-band light curves reliably recover the mass ratios for contact binary systems determined independently using spectroscopy as long as the ratio of stellar temperatures is constrained to near unity (our T 2 ≈T 1 model). Temperature ratio T 2 /T 1 and mass ratio q can be degenerate in ways that yield incorrect mass ratios if T 2 /T 1 is unconstrained. Light curve modeling combined with phase-and velocity-resolved spectroscopy-in particular, the broadening function analysis-provides strong constraints on system parameters for the ten near-(contact)systems having both datasets. Adopted practices include two rounds of MCMC model fitting, the first to locate the global minimum within the multi-modal parameter space and the second to characterize the shape and extent of the minimum after marginalizing over the (unknown and unmodeled) parameters by re-normalizing the reduced chi-squared of the best-fitting model to one. After experimenting with a small suite of competing models, we adopt either a five-free-parameter contact model or a six-free-parameter detached model. Even when additional model components such as a spot are required (KIC08913061, KIC09164694, KIC09345838) the mass ratio is recovered with good fidelity. The MCMC 1σ uncertainties are generally small for log q and cos i while f and l3 are still constrained but with less precision. 2. The PHOEBE model line profiles agree reasonably well with the broadening functions determined from the R≈4000 optical spectra obtained near quadrature phases. The differences between the model line profiles and the broadening functions provide measures of the third-light fractions that are usually in good agreement with those independently obtained from the best-fitting models. The broadening functions in four of the ten spectroscopically studied systems (KIC04999357, KIC06844489, KIC10292413, and KIC11097678) show a deficit relative to the model at velocities associated with the secondary (less luminous) component. Such a deficit could be caused by cool spot or line-emitting material on the secondary. We favor the former possibility because the deficit persists when portions of the spectrum such as Hα are excluded from the BF analysis. These BF deficits in extreme-q may also be explained in terms of accretion flows between components, as in AW UMa (Rucinski 2015) . The characteristics of persistent surface inhomogeneities on one or both components could be further constrained using time-resolved multi-band photometry and/or spectroscopy. 3. At least eight of the ten spectroscopically studied contact binaries show evidence for a third component, either in the O − C residuals or the broadening function, consistent with prior evidence that most contact binaries exist in triple systems (Pribulla & Rucinski 2006 ). Horvat et al. (2019) implementation of irradiation effects yields mixed results on the ten pilot study systems. Sometimes including irradiation in the models leads to a significantly lower RMS (KIC04999357), and sometimes the RMS is the same or larger (most of the pilot sample). Often the impact of spots and/or intrinsic variability swamps the effects of irradiation (KIC08913061, KIC09164694, KIC09345838). Generally, inclusion of irradiation effects changes the best-fitting parameters very little. Among the probable contact systems T 2 /T 1 becomes slightly less extreme and mass ratios become slightly less extreme when irradiation effects are included in the models. Given the ubiquity of photometric variability arising from spots and other stellar activity, these effects may mask or be degenerate with those caused by temperature differences and irradiation. We conclude that implementing irradiation effects makes a minimal or negligible difference in the retrieval of parameters for the majority of contact binary systems. Detached systems having large temperature differences may show otherwise. 5. For detached binaries evincing large ellipsoidal modulations like KIC09953894, we show that it is possible to retrieve approximate system parameters, including q, without spectroscopic knowledge of the mass ratio. The allowed range of parameters is larger, however, than when q is known a priori. 6. Contact binary models of 783 Kepler systems with orbital periods P <2 d yield 178 (23%) that are best reproduced using the T 2 ≈T 1 contact configuration. Cos i and log q are generally well constrained (median 1σ ranges 0.045 and 0.10, respectively) while fillout factors and third-light fractions are less certain (1σ ranges 0.15 and 0.08, respectively) and often have some degree of degeneracy with other parameters. Constraints are less strong among the 114 (15%) of systems deemed detached binaries. Over half (62%) of modeled systems were classified as ambiguous geometries; most of these have best-fitting models with large RMS resulting from asymmetric light curves or best-fitting models requiring unequal temperature components. From single-band light curves alone it is not possible to distinguish a low-fillout-factor contact system from a detached configuration. Especially for systems having P <0.5 d the RMS for the best-fitting contact and detached models are often similar. At these short periods the best-fitting solutions from contact and detached models tend to yield very similar parameters, but at longer periods the solutions can differ considerably. 7. The distribution of third-light fractions presented in Figure 60 may be the most robust characterization of this parameter yet achieved, given the high precision of the Kepler photometry coupled with the unbiased and complete nature of the Kepler contact binary sample. We find that 77% of contact binaries have third light fractions l3>0.15, meaning that bright tertiaries are common among contact binary systems and that the true fraction of triple systems among this population is even larger. 8. The vast majority of systems at long periods P >0.5 d have either detached or ambiguous classifications, indicating that true contact binaries become rare at these longer periods. The Kirk et al. (2016) light curve morphology parameter serves as a rough indicator of the geometrical configuration of the close binaries studied (0.700.9 and long periods, affirming the consistency of these independent approaches to characterizing physical geometry through analysis of the light curves. 9. Among probable contact systems, the distribution of fillout factors covers the full range from 0-1 but shows a small excess near 1, indicating a measurable subpopulation with Roche lobes near the maximum permitted volumes. 10. Mass ratios q a >0.8 are rare among the probable contact binary subset and the ambiguous subset, providing evidence for rapid evolution away from mass ratios near unity, as shown in Figure 63 (Section 5.2) from the new set of model predictions of Molnar (2022) . The distribution of systems in q and P is consistent with a scenario whereby conservative evolution drives systems toward longer periods and more extreme mass ratios, more rapidly for systems that begin with mass ratios near unity than for systems that begin with small mass ratios. Mass ratios q a <0.1 are rare among the whole sample and appear to be absent among contact binaries, in agreement with binary evolution models which predict a critical lower limit as a function of period, below which the onset of a Darwin instability promotes a rapid merger. Pre-merger systems should be found at this lower boundary with periods P >0.75 d and q=0.045-0.15. These conclusions affirm the power of modern binary modeling codes combined with MCMC techniques to retrieve system parameters for large populations of close binaries from single-band light curves. Although the computational effort required is considerable ( 80,000 CPU days for this project involving 783 systems) suitable resources are becoming common even outside of major supercomputer centers. The addition of photometric data at a second bandpass would provide further constraints on third-light parameters and reduce degeneracies owing to T 2 /T 1 , thereby narrowing the dispersion of posterior parameter distributions. Multi-band high-cadence long-baseline datasets are imminent, if not already available. When coupled with planned time-domain photometric and spectroscopic surveys, the methods demonstrated here are a potent tool when applied in a statistical manner to reveal the evolutionary pathways and role of distinct physical processes (Lidov-Kozai cycles, tidal friction, mass transfer, magnetic cycles, thermal relaxation oscillations, common envelope evolution, etc.) among large populations of close binaries. Note-(1)-Kepler Input Catalog identifier, (2)-mean Kepler band magnitude, (3)-stellar effective temperature from the Kepler Input Catalog, (4)-light curve morphology parameter from Kirk et al. (2016) , (5)-orbital period determined as described in this work, (6)reference time of superior conjunction, as determined in this work, e Best-fitting model involves a hot spot on the primary. f A visual companion at 2 separations and ∆G=0.5 mag may contribute light within the Kepler pixel. g The secondary is very faint and its radial velocity is significantly uncertain. Table 3 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. A (1) in the final column designates systems having superior model fits using detached configuration parameters summarized in Table 5 . Table 4 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. Table 5 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. R1/R1max and R2/R2max are approximate values predicated on the assumption of M =1 M for the more massive component. R1max and R2max scale weakly with the adopted M1. R1 is the radius of a spherical star having equivalent surface area to the tidally distorted star. R1max is the equivalent radius of a tidally distorted star of maximum possible size without overflow. Table 6 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. R1/R1max and R2/R2max are approximate values predicated on the assumption of M =1 M for the more massive component. R1max and R2max scale weakly with the adopted M1. R1 is the radius of a spherical star having equivalent surface area to the tidally distorted star. R1max is the equivalent radius of a tidally distorted star of maximum possible size without overflow. Modelling of Stellar Atmospheres The Observatory Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series American Astronomical Society Meeting Abstracts, American Astronomical Society Meeting Abstracts American Astronomical Society Meeting Abstracts #233 Modeling and Analysis of Eclipsing Binary Stars; The theory and design principles of PHOEBE Astronomical Society of the Pacific Conference Series The Zwicky transient facility observing system Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series We thank undergraduate students Elle Buser and Alexander Schultz for helping obtain data in support of this program. We are grateful to Adam Myers for consultation on the use of MCMC methods and the emcee code. We thank Natalie Batalha for communications regarding the selection of the Kepler target sample. We thank Andre Prśa and the PHEOBE 2 team for making a powerful tool available to the community. This work made use of the Advanced Research Computing Center (2018) Teton Computing Environment, Intel x86 64 cluster, University of Wyoming, Laramie, WY https://doi.org/10.15786/M2FY47. This paper was made possible, in part, by the extended isolation enforced by the 2020-2021 SARS-CoV2 pandemic. It is dedicated to Eugene M. Hermansen (U.S. Army, 1951Army, -1953, who became one of the casualties during the 2020-2022 pandemic. This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research was supported by the National Science Foundation through grant AST-1716718. Software: IRAF (Tody 1986) , PHOEBE , emcee (Foreman-Mackey et al. 2013)