Constraining the Circumbinary Disk Tilt in the KH 15D system
MMNRAS , 1–15 (2020) Preprint 1 October 2020 Compiled using MNRAS L A TEX style file v3.0
Constraining the Circumbinary Disk Tilt in the KH 15Dsystem
Michael Poon, , (cid:63) J. J. Zanazzi, † and Wei Zhu, ‡ Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, Ontario, M5S 1A7, Canada Department of Astronomy and Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
KH 15D is a system which consists of a young, eccentric, eclipsing binary, and acircumbinary disk which obscures the binary as the disk precesses. We develop a self-consistent model that provides a reasonable fit to the photometric variability that wasobserved in the KH 15D system over the past more than 60 years. Our model suggeststhat the circumbinary disk has an inner edge r in (cid:46) , an outer edge r out ∼ a few au ,and that the disk is misaligned relative to the stellar binary by ∼ ∼
10 degrees and ∼
15 degrees, respectively. We also provide constraints onother properties of the disk, such as the precession period and surface density profile.Our work demonstrates the power of photometric data in constraining the physicalproperties of planet-forming circumbinary disks.
Key words: protoplanetary discs – planets and satellites: formation – stars: individ-ual (KH 15D) – binaries: spectroscopic – techniques: photometric
Our understanding of tilts within planet-forming circumbi-nary systems has undergone drastic changes within the pastdecade. Originally, the basic picture was quite simple: a cir-cumbinary disk should always be observed to be alignedwith the orbital plane of the binary. Even though simula-tions of turbulent molecular clouds found circumbinary disksfrequently formed misaligned with the orbital plane of thebinary (e.g. Bate et al. 2002; Bate 2012, 2018), viscous disk-warping torques were showed to damp the disk-binary incli-nation over timescales much shorter than typical protoplan-etary disk lifetimes (Foucart & Lai 2013, 2014). Most in-clination constraints on protoplanetary (e.g. Andrews et al.2010; Rosenfeld et al. 2012; Czekala et al. 2015, 2016; Ru´ız-Rodr´ıguez et al. 2019) and debris (e.g. Kennedy et al. 2012b)disks confirm this basic picture, finding alignment of the diskwith the orbital plane of the binary to within a few degrees.However, after the detection of a few highly inclinedcircumbianry disks (e.g. Kennedy et al. 2012a; Marino et al.2015; Brinch et al. 2016; Czekala et al. 2017), it became clearnot all circumbinary disks align rapidly. Motivated by these (cid:63) [email protected] (MP) † [email protected] (JJZ) ‡ [email protected] (WZ) detections of highly-inclined disks, the theoretical commu-nity found that when a circumbinary disk orbits an eccentric binary, the disk-binary inclination can grow under certaincircumstances, evolving eventually to ◦ (polar alignment;Aly et al. 2015; Martin & Lubow 2017; Zanazzi & Lai 2018).Additional inclined circumbinary disks orbiting eccentric bi-naries were discovered soon thereafter, such as HD 98800(Kennedy et al. 2019), and AB Aurigae (Poblete et al. 2020).Recently, Czekala et al. (2019) showed circumbinary diskshad larger inclinations when orbiting binaries with highereccentricities, further supporting the operation of this mech-anism in circumbinary disk systems.The alignment process itself was also shown to be non-trivial, with the disk itself occasionally breaking in the pro-cess. Early on, the disk was expected to remain nearly flat,due to the resonant propagation of bending waves acrossthe disk (Papaloizou & Lin 1995; Lubow & Ogilvie 2000).Later hydrodynamical simulations found that the disk un-der some circumstances may break, with different disk an-nuli becoming highly misaligned with one another, due tostrong differential nodal precession induced by the torquefrom the binary (e.g. Nixon et al. 2013; Facchini et al. 2013).Numerous broken protoplanetary disks orbiting two binarystars have subsequently been found, including HD 142527(Marino et al. 2015; Price et al. 2018) and GW Ori (Bi et al.2020; Kraus et al. 2020). © a r X i v : . [ a s t r o - ph . E P ] S e p Poon, Zanazzi, & Zhu
The inclinations of detected circumbinary planets,in contrast, remain broadly consistent with formation innearly-aligned circumbinary disks. After the detection of afew dozen circumbinary planets (see Welsh & Orosz 2018;Doyle & Deeg 2018 for recent reviews), the inclinationswithin the circumbinary planet population are consistentwith alignment to within ∼ ◦ (Martin & Triaud 2014; Arm-strong et al. 2014; Li et al. 2016). However, highly-misalignedcircumbinary planets are physically allowed, because the in-clined orbit has been shown to be long-term stable once aplanet forms in the polar-aligned circumbinary disk. (Doolin& Blundell 2011; Giuppone & Cuello 2019; Chen et al. 2019,2020). New detection methods may detect polar-aligned cir-cumbinary planets in the future (Zhang & Fabrycky 2019).While a large number of systems now have constraintson mutual inclinations between the disk and the binary or-bital planes, there remain few constraints on twists andwarps within the circumbinary disk itself. This is becausethe methods used to constrain disk inclinations are not sen-sitive to the small misalignments within the disk. Gaseousprotoplanetary disk inclinations are constrained via the or-bital motion of the disk gas through the Doppler shift ofemission lines (e.g. Facchini et al. 2018; Price et al. 2018).Debris disk inclinations are constrained by the orientationof the disk implied by its continuum emission (e.g. Kennedyet al. 2012a,b).A rare example of a circumbinary disk system wherephotometric constraints exist is Kearns-Herbst 15D (KH15D) (Kearns & Herbst 1998). KH 15D is a system witha highly-unusual light-curve, which exhibited dips by up to5 orders of magnitude. The morphology of the dipping be-havior changed over decade-long timescales, but displayedperiodicity over short 48 day timescales. The complex light-curve of this system is generally believed to be due to ancircumbinary disk and an eclipsing binary, with some of thedips caused by the optically thick, precessing disk slowly ob-scuring the orbital plane of the stellar binary (Winn et al.2004, 2006; Chiang & Murray-Clay 2004; Capelo et al. 2012).But although much work has gone into understanding KH15D, no work has attempted to provide quantitative con-straints on the properties of the warped disk based on thephotometric data.In this work, we combine the spectroscopic and photo-metric data to constrain the properties of the circumbinarydisk KH 15D. With recent data up to 2018 from Aronowet al. (2018) and Garc´ıa Soto et al. (2020), we improve theWinn et al. (2006) model to fit all photometric data since1955. Our results are particularly exciting, as our fit nearlyencompasses the full transit of the circumbinary disk of KH15D. Section 2 extends the Winn et al. (2006) model to fitthe light-curve of the system, over the more than 60 year du-ration the system was observed. Section 3 develops a dynam-ical model to constrain the circumbinary disk properties im-plied by the photometric constraints. Section 4 discusses thetheoretical implications of our work, improvements whichcan be made to our model, and our model predictions whichcan be tested with future observations. Section 5 summarizesthe conclusions of our work. We use radial velocity observations to constrain the orbitof the stellar binary and photometric data to constrain thegeometry of the optically thick, precessing disk. We beginwith a brief description of the radial velocity and photomet-ric data used in this work.We use the Radial Velocity (RV) measurements gath-ered by Hamilton et al. (2003) and Johnson et al. (2004),and because one of the stars is eclipsed by the disk whenthe RV measurements are taken, this is effectively a single-lined spectroscopic binary. As in Winn et al. (2006), weonly use RV measurements gathered when the system fluxis 90% or greater than its mean out-of-eclipse flux, becausethe Rossiter-McLaughlin effect (Rossiter 1924; McLaughlin1924) leads to systematic errors as the stellar companion isocculted by the disk. This gives 12 RV measurements to aidin constraining the orbit of the binary.For the photometry, we use the tabulated data fromWinn et al. (2006), Aronow et al. (2018), and Garc´ıa Sotoet al. (2020). Details about the photometric observations canbe found in these references, and we only provide a brief sum-mary here. These catalogues include data from photographicplates from the 1950s to 1985 (Johnson & Winn 2004; Maffeiet al. 2005), as well as observations using Charge-CoupledDevices (CCDs) since 1995 (e.g. Hamilton et al. 2005; Capeloet al. 2012; Aronow et al. 2018; Garc´ıa Soto et al. 2020). Noobservations were known to be taken between 1985-1994.All photometric observations have been transformed into thestandard Cousin I -band measurements. We bin the originaldata-set (6241 points) into 2813 data points in order to re-duce the amount of time needed for the photometric modelcomputation. The uncertainties of all photometric measure-ments are re-scaled up by a factor of two to produce mean-ingful statistical metrics (i.e., reduced χ = ). So far three models have been proposed to explain the pho-tometric variation in the KH 15D system. The phenomeno-logical model by Winn et al. (2004, 2006) approximates theleading edge of the disk as an infinitely long and opticallythick screen, which eclipses the two stars as the screen movesacross the orbit of the binary. Motivated by dynamics, Chi-ang & Murray-Clay (2004) treat the KH 15D disk as awarped disk with finite optical depth, and model the pho-tometric variations by a disk precessing into and out of theline-of-sight of the observer. Silvia & Agol (2008) developedtheir model based on the model of Winn et al. (2006), butintroduced more disk-related physics, such as the finite op-tical depth, curvature near the edge, and forward-scatteringof starlight from the dust in the disk (which was parameter-ized as “halos” in the Winn et al. 2006 model). We chooseto build our model based on Winn et al. (2006), because itallows us to remain agnostic about the detailed physics ofthe disk itself, while still accurately fitting the light-curveof KH 15D. We review the Winn et al. (2006) model withinthis subsection.Figure 1 illustrates the physical motivation behind the
MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Circumbinary Disk
Visualization
Star B
Y X ω Y(t) θ(t)
Star A zoomed-in
Figure 1.
The purple ellipses displays a potential shape for the circumbinary disk of KH 15D. No assumption about the inclinationbetween the innermost edge (red dashed curve), outermost edge (yellow dashed curve) or binary plane (centre yellow ellipses) are madewhile modelling the eclipses of KH 15D. We note that our model also allows for the leading edge to be the outermost truncation radius ofthe disk, with the tailing edge the innermost truncation radius of the disk. The zoomed-in inset diagram displays how the circumbinarydisk geometry eclipses the binary of KH 15D. The inner or outer truncation radius of the disk slowly covers the orbital plane of thebinary, as the disk precesses around the orbital angular momentum axis of the binary. We approximate the inner and outer disk edgesas straight edges as the binary is eclipsed.
Winn et al. (2006) model. A single “leading” edge (reddashed line) slowly advances over the orbital plane of thebinary, which approximates the inner or outer truncationradius of the disk slowly covering both stars as the disk isprecessing around the binary. For data taken before 2005,it is reasonable to neglect the outer disk truncation radius(the “tailing” edge, as marked by the yellow dashed curve).To model how quickly the leading edge advanced across theorbit of the binary, Winn et al. (2006) used the latest datewhen star B was still visible ( t ), and the latest date whenthe orbit of star A was visible ( t , see Fig. 2 left panel), asfree parameters in his photometric model. In addition, notonly was the angle the leading edge made with the X -axisof the observer allowed to vary, but it was also allowed tochange at a constant rate, controlled by two free parame-ters θ L ( t ) and (cid:219) θ L . The light from the binary was modelledwith 7 parameters, with the luminosity from star A (B) de-noted by L A ( L B ), the background light when the disk fullyeclipses the binary by L , with the parameters { (cid:15) , (cid:15) , ξ , ξ } parameterizing the light emitted by halos surrounding starsA and B. Specifically, the 1D brightness distribution fromstar i = A , B was taken to be B i ( v ) = ( (cid:15) / ξ ) exp [( v + )/ ξ ] v ≤ − ( (cid:15) / ξ ) + B (cid:63) i ( v ) − < v < ( (cid:15) / ξ ) exp [−( v − )/ ξ ] v ≥ , (1)where B (cid:63) i is the 1D brightness distribution of star i , assum-ing a linear limb-darkening model: B (cid:63) i ( v ) = I i (cid:112) − v (cid:104) − u (cid:16) − π (cid:112) − v (cid:17)(cid:105) . (2)Here, u = . is the limb-darkening coefficient for bothstars, and I i is the reference intensity of star i . Letting y L , i bethe distance of the lead edge from star i , and v L , i = y L , i / R i , then the flux from star i is F L , i = ∫ ∞ v L , i B i ( v ) d v , (3)with the total flux F = F L , A + F L , B . Physically, each “halo”parameterizes forward-scattering of starlight by dust in thedisk (Winn et al. 2006; Silvia & Agol 2008). The massand radius for star A were taken to be M A = . (cid:12) and R A = . (cid:12) , while the ratios between the masses and radiiof the two stars are M B / M A = . and R B / R A = . , re-spectively. The orbit of the binary is described by standardorbital parameters used to model RV data, with an orbitalperiod P , eccentricity e , inclination I , longitude of pericenter ω , time of pericenter passage T p , and line-of-sight velocity γ (see e.g. Fulton et al. 2018 for details). The Cartesian co-ordinate system in the sky-projected reference plane of theobserver ( X , Y ) is chosen so the X -axis lies along the line ofnodes (so Ω = ).The best fit model parameters for the KH 15D systemwas then calculated by minimizing (Winn et al. 2004, 2006) χ = N F (cid:213) j = (cid:18) F j − F O , j σ F , j (cid:19) + λ N V (cid:213) j = (cid:18) V j − V O , j σ V , j (cid:19) ≡ χ + λ χ , (4)where χ is the χ of the photometry model alone, χ isthe χ metric of the modelled orbit of the binary in relationto the RV data (see Fulton et al. 2018 for details), and for aquantity X , X j denotes the model prediction at point j , X O , j denotes the observed value of X at j , while σ X , j denotes theuncertainty of X O , j at j . The parameter λ = was chosen toincrease the importance of the RV model relative to that forthe photometry, because the model constraining the orbitof the binary (a Keplerian orbit) is much more certain than MNRAS000
Winn et al. (2006) model. A single “leading” edge (reddashed line) slowly advances over the orbital plane of thebinary, which approximates the inner or outer truncationradius of the disk slowly covering both stars as the disk isprecessing around the binary. For data taken before 2005,it is reasonable to neglect the outer disk truncation radius(the “tailing” edge, as marked by the yellow dashed curve).To model how quickly the leading edge advanced across theorbit of the binary, Winn et al. (2006) used the latest datewhen star B was still visible ( t ), and the latest date whenthe orbit of star A was visible ( t , see Fig. 2 left panel), asfree parameters in his photometric model. In addition, notonly was the angle the leading edge made with the X -axisof the observer allowed to vary, but it was also allowed tochange at a constant rate, controlled by two free parame-ters θ L ( t ) and (cid:219) θ L . The light from the binary was modelledwith 7 parameters, with the luminosity from star A (B) de-noted by L A ( L B ), the background light when the disk fullyeclipses the binary by L , with the parameters { (cid:15) , (cid:15) , ξ , ξ } parameterizing the light emitted by halos surrounding starsA and B. Specifically, the 1D brightness distribution fromstar i = A , B was taken to be B i ( v ) = ( (cid:15) / ξ ) exp [( v + )/ ξ ] v ≤ − ( (cid:15) / ξ ) + B (cid:63) i ( v ) − < v < ( (cid:15) / ξ ) exp [−( v − )/ ξ ] v ≥ , (1)where B (cid:63) i is the 1D brightness distribution of star i , assum-ing a linear limb-darkening model: B (cid:63) i ( v ) = I i (cid:112) − v (cid:104) − u (cid:16) − π (cid:112) − v (cid:17)(cid:105) . (2)Here, u = . is the limb-darkening coefficient for bothstars, and I i is the reference intensity of star i . Letting y L , i bethe distance of the lead edge from star i , and v L , i = y L , i / R i , then the flux from star i is F L , i = ∫ ∞ v L , i B i ( v ) d v , (3)with the total flux F = F L , A + F L , B . Physically, each “halo”parameterizes forward-scattering of starlight by dust in thedisk (Winn et al. 2006; Silvia & Agol 2008). The massand radius for star A were taken to be M A = . (cid:12) and R A = . (cid:12) , while the ratios between the masses and radiiof the two stars are M B / M A = . and R B / R A = . , re-spectively. The orbit of the binary is described by standardorbital parameters used to model RV data, with an orbitalperiod P , eccentricity e , inclination I , longitude of pericenter ω , time of pericenter passage T p , and line-of-sight velocity γ (see e.g. Fulton et al. 2018 for details). The Cartesian co-ordinate system in the sky-projected reference plane of theobserver ( X , Y ) is chosen so the X -axis lies along the line ofnodes (so Ω = ).The best fit model parameters for the KH 15D systemwas then calculated by minimizing (Winn et al. 2004, 2006) χ = N F (cid:213) j = (cid:18) F j − F O , j σ F , j (cid:19) + λ N V (cid:213) j = (cid:18) V j − V O , j σ V , j (cid:19) ≡ χ + λ χ , (4)where χ is the χ of the photometry model alone, χ isthe χ metric of the modelled orbit of the binary in relationto the RV data (see Fulton et al. 2018 for details), and for aquantity X , X j denotes the model prediction at point j , X O , j denotes the observed value of X at j , while σ X , j denotes theuncertainty of X O , j at j . The parameter λ = was chosen toincrease the importance of the RV model relative to that forthe photometry, because the model constraining the orbitof the binary (a Keplerian orbit) is much more certain than MNRAS000 , 1–15 (2020)
Poon, Zanazzi, & Zhu the model describing the light-curve of the binary (eclipsedby a precessing disk).Because the screen advances in the positive vertical di-rection at a constant rate, an equivalent way of parameter-izing the ascent of the screen are through where the screenintersects the Y -axis at the orbital contact time t , whichwe will denote by Y L ( t ) (see Figs. 1 & 2), and the rate ofchange in the Y -direction, (cid:219) Y L . Winn et al. (2006) choose t and t because of its tighter connections with observations( t , t denote changes in the light-curve of KH 15D). Whenextending the Winn et al. (2006) model, we will also primar-ily refer to orbital contact times to parameterize the advanceof the screen across the orbit of the binary (Fig. 2), but alsofrequently refer to Y L and (cid:219) Y L as well. Our photometric model builds off of Winn et al. (2006), andseeks to fit the light-curve of KH 15D from 1955-2018 withminimal modifications. After 2012, the tailing edge startedto uncover star B, due to the other (inner or outer) trunca-tion radius of the disk of KH 15D precessing over the binaryorbit with respect to the line-of-sight of the observer (seeFig. 1). The simplest extension is to include an additionaltailing edge in the modelling (denoted by subscript T ), whichlags in position behind the leading edge (with subscript L ),which intersects the Y -axis at a location Y T ( t ) (see Fig. 2).This tailing edge also introduces 5 new orbital contact timesas the edge crosses the orbit of the binary: t , t , t , t , t (seeFig. 2 for illustration). Assuming θ L = θ T and (cid:219) Y L = (cid:219) Y T , theprevious 1-edge semi-infinite sheet becomes a 2-edge thinrectangular sheet of constant width, which is infinite alongits length. Garc´ıa Soto et al. (2020) assumed this for theirlight-curve model, and neatly fit CCD photometry from 1995and onwards. However, because this fit does not match thelight-curve data prior to 1995 (fit not shown here or in Gar-c´ıa Soto et al. 2020), further modifications are needed to theWinn et al. (2006) and Garc´ıa Soto et al. (2020) model.Through much experimentation, we found the followingset of additions to the Winn et al. (2006) model that let us fitthe 60+ year light-curve. The connection of these additionswith a warped disk driven into precession by an eccentricbinary will be made clear in the following section. • We let the leading and tailing edges have different an-gles ( θ L [ t ] (cid:44) θ T [ t ] , see Fig. 2). • We let each edge linearly evolve in time independently( (cid:219) θ L [ t ] (cid:44) (cid:219) θ T [ t ] ). • Parameterize (cid:219) θ L by two constant, piecewise rates intime: (cid:219) θ L ( t ) = (cid:219) θ L when t < t , and (cid:219) θ L ( t ) = (cid:219) θ L when t > t .We keep (cid:219) θ T ( t ) = constant as a single parameter. Because wemake the leading edge symmetric about t , we fit for thetimes { t , t } in our MCMC model to constrain Y L ( t ) and (cid:219) Y L , rather than { t , t } as in Winn et al. (2006). • Let the width of the screen to change over time ( (cid:219) Y L (cid:44) (cid:219) Y T ), but keep both rates (cid:219) Y L and (cid:219) Y T constant with time. • Prescribe the rate of ascent of the tailing edge in rela-tion to the rate of ascent of the leading edge. Specifically,we take (cid:219) Y T = α (cid:219) Y L for α = . , . , . , . , . , . . We alsoexperimented with letting (cid:219) Y T be a free parameter (fitting forthe contact times { t , t } ), and found these fits gave (cid:219) Y T ≈ (cid:219) Y L ,but the MCMC did not always converge. We choose this pa- Table 1.
Definitions of Model Parameters.Free Parameter Description P Orbital period e Orbital eccentricity I Inclination of orbital plane ω Argument of pericenter T p Time of periapsis passage L B / L A Luminosity of star B relative to star A (cid:15) Fractional flux of stellar halo (cid:15) Fractional flux of stellar halo ξ Exponential scale factor of stellar halo ξ Exponential scale factor of stellar halo t Third orbital contact time t Fifth orbital contact time t Sixth orbital contact time θ L ( t ) Angle between x-axis and leading edge at t = t θ T ( t ) Angle between x-axis and tailing edge at t = t (cid:219) θ L Rotation rate of leading edge when t < t (cid:219) θ L Rotation rate of leading edge when t > t (cid:219) θ T Rotation rate of tailing edge In the direction the leading edge approaches the star In the direction the leading edge travels beyond the star Defined in Figure 2 rameterization to make sure the other model parameters arewell-determined.We further simplify the Winn model by analytically solvingfor L A and γ with respect to the rest of the parameters,since they are constant shifts to the photometric and radialvelocity models, respectively. This reduces the number offree parameters by 2. For reference, we display each modelparameter and its definition in Table 1.To calculate the flux from the KH 15D system, we sim-ply add the flux from stellar light emitted exterior to thetailing edge, to that emitted exterior to the leading edge. Inmore detail, letting y T , i be the distance of the tailing edgefrom star i , with v T , i = y T , i / R i , star i emits the flux F T , i = ∫ v T , i −∞ B i ( v ) d v (5)exterior to the tailing edge, giving the total flux F = F L , A + F L , B + F T , A + F T , B . We neglect the intersection between theleading and tailing edges in the flux calculation, because thisintersection occurs far from the orbit of the binary.For our radial velocity model, we follow equations (2)and (3) from Section 2.1 of Fulton et al. (2018). To opti-mize the model parameters, we use a Python-implementedMarkov chain Monte Carlo (MCMC) package emcee byForeman-Mackey et al. (2013). We use the same χ statisticas in equation (4).Preliminary tests find the background light in the KH15D system to be L ≈ , so we remove L from our modelparameters. This is expected if L is from forward scatteringof the stars’ light around the tailing edge of the disk (Silvia& Agol 2008), rather than the finite optical depth of the diskitself (Chiang & Murray-Clay 2004), because forward scat-tering of stellar light around both screen edges is included inour model. Our final model has 18 parameters (see Table 1),which we run for 20,000 steps with 36 walkers. Running thefinal model for each α , we come to the following results: eachMCMC converged except for α = . , with only α = . , . MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Figure 2.
Left Panel : Definitions of orbital contact times. The leading (tailing) edge has contact times t to t ( t to t ). All contacttimes denote when the leading or tailing edge lies tangent to the orbit of either star A or B, with the exception of t and t , whichdenote when the leading or tailing edges intersect the center of mass of the binary. Right Panel : The definitions of quantities related toour model of the disk eclipsing binary of KH 15D, which we model as an opaque screen bounded by two infinitely-long, straight edgeson both sides. The leading (tailing) edge is parameterized by its intersection with the Y -axis, Y L ( Y T ), and the angle between the edgeand the X -axis, θ L ( θ T ). Our model allows for Y L , θ L , Y T , and θ T to evolve (linearly) with time. We emphasize our model makes noassumption on the underlying geometry of the disk eclipsing the binary of KH 15D. producing reasonable-looking light-curves. Model parame-ters for α = . , . , . , . are reported in Appendix A1. Wehighlight α = . as the best fit with parameters in Table 2,and display corner plots of the posteriors in Appendix B1.Our fit for the entire light-curve of KH 15D is displayedin Figure 3. Our model does a good job in describing boththe maximum and minimum fluxes from KH 15D, whichchange with time. As expected, the orbital contact times t i denote when the light-curve of the system changes its mor-phology. The gradual change of maximum/minimum fluxaround t i values is due to the halos around each star: forpoint-source stars eclipsed by a razor-thin opaque edge, thephotometric model predicts almost discontinuous changes inlight-curve morphologies around t i values. Figures 4-5 show the photometric model and data foldedover the binary orbital period, which is comparable to Fig-ure 12 in Winn et al. (2006). Again, we see our model doesa good job at modelling changes in the light-curve of KH15D, with orbital contact times (see Table 2 for values) de-lineating morphology changes as one or both stars becomeseclipsed or revealed by an edge. Examining the data from2013-2014 and 2015-2016, the large scatter makes it seemunlikely that any (simple) model could provide an accuratefit to the observed light-curve. In addition, we do not removeany outliers (compare 1998-1999 panel in Fig. 4 to Fig. 12in Winn et al. 2006). An interesting feature occurs around2010, when the egress is poorly fit (for all values of α ). Thiscould be related to the clumpiness/transparency near theedges of the disk as discussed in Garc´ıa Soto et al. (2020),where the assumption of sharp edges breaks down. The previous section showed that in order for the Winn et al.(2006) model to fit the entire more than 60 year light-curveof KH 15D, a number of modifications to this original modelmust be made. In this section, we illustrate how these mod-ifications are motivated by the dynamics of a warped disk,driven into precession around an eccentric binary. In doingso, we will show that the disk of KH 15D orientation, warp,radial extent, and even surface density profile may be con-strained by photometry alone.
MNRAS000
MNRAS000 , 1–15 (2020)
Poon, Zanazzi, & Zhu I [ m a g ] F l ux t t t t t t t Figure 3.
Light-curve of KH 15D from 1965 to present, displaying the complex change in variability seen with time. The observed lightcurve in I -band magnitude is shown in the upper panel, and the light curve after the normalization to the flux of star A is shown inthe lower panel. Blue points are photometry from Aronow et al. (2018) and Garc´ıa Soto et al. (2020), while the thin black line displaysour photometric model fit (see Table 2 for parameter values). Vertical cyan lines denote the orbital contact times t i indicated, wherethe leading or tailing edge of the screen (e.g. circumbinary disk, see Fig. 1) hits a different portion of the binary orbit (see Fig. 2 fordefinitions). Our model does well in reproducing the KH 15D light-curve variability over the length of time the system is observed. For the disk around KH 15D to coherently precess over itslifetime, internal forces within the disk must keep neighbor-ing disk annuli nearly aligned with one another, otherwisedifferential nodal precession from the gravitational influenceof the binary will disrupt and “break” the disk (e.g. Larwood& Papaloizou 1997; Facchini et al. 2013; Nixon et al. 2013;Martin & Lubow 2018). When these internal torques aremuch stronger than the external torque on the disk from thebinary, the disk behaves as a rigid body, coherently precess-ing about the orbital angular momentum axis of the binary(e.g. Martin & Lubow 2017; Smallwood et al. 2019; Moodyet al. 2019). To model the dynamical evolution of the disk,we will assume the disk behaves approximately like a rigidplate, treating the disk as a secondary whose mass is dis-tributed between radii r L and r T .However, before we introduce our model for an extendeddisk, we discuss the dynamics of a test particle on a circularorbit (which we will refer to as a ring), driven into precessionby the torque from the binary. Many authors have shown theorbital angular momentum unit vector of the ring ˆ l r is driveninto precession and nutation about either the orbital angu-lar momentum unit vector of the binary ˆ l , or eccentricityvector of the binary e (vector in pericenter direction with magnitude e ). The dynamical evolution of the ring dependssensitively on the initial orientation of ˆ l r with respect to ˆ l and e , as well as the magnitude of the eccentricity of thebinary e . To calculate the evolution of ˆ l r about ˆ l and e , weadopt the formalism of Farago & Laskar (2010), who calcu-lated the secular evolution of a ring about a massive binarywith an eccentric orbit, after expanding the Hamiltonian ofthe binary to leading order in r / a (where r is the semi-majoraxis of the test particle), and averaging over the mean mo-tions of the test particle and the binary. It was found thecharacteristic precession and nutation frequency of ˆ l r aboutthe binary was given by (denoted by α in Farago & Laskar2010) ν = µ M t (cid:18) GM t a (cid:19) / (cid:16) ar (cid:17) / , (6)where M t = M A + M B is the total mass of the binary, while µ = M A M B / M t is the reduced mass of the binary.After calculating the evolution of a (circular) test par-ticle ˆ l r vector about ˆ l and e using Farago & Laskar (2010),we then translate the evolution of ˆ l r into the inclination ofthe test particle I r and longitude of ascending node Ω r , inthe frame where ˆ z = ˆ l and the line of nodes points in thedirection of e . Because the orientation of the binary orbit inthe reference frame of a distant observer is described by the MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Figure 4.
Data (blue points) and fitted model (black lines) displayed in Figure 3, folded over the binary orbital period, prior to theyear 2005. The timespan over which the data and model are folded over is displayed in each figure. Our model reproduces the changingmorphology of the light-curve of KH 15D well. orbital elements { a , e , ω, I , Ω } , the position of the ring in theframe of the observer is (cid:169)(cid:173)(cid:171) x y z (cid:170)(cid:174)(cid:172) r , obs = R Z ( Ω ) R X ( I ) R Z ( ω ) R Z ( Ω r ) R X ( I r ) (cid:169)(cid:173)(cid:171) x r y r (cid:170)(cid:174)(cid:172) , (7)where ( x r , y r ) = r ( cos ϕ, sin ϕ ) parameterizes the ( X , Y ) coor-dinates of the ring in the frame where ˆ z = ˆ l , and R X [ β ] ( R Z [ γ ] ) denote rotations along the X ( Z )-axis by angles β ( γ ). As in Winn et al. (2006), we choose the reference planeof the observer so the X -axis points along the binary line-of-nodes (so Ω = ). Also, because our MCMC model highly favors a nearly edge-on orbit (Table 2), we assume I (cid:39) ◦ for simplicity for the rest of this section. All other orbitalparameters are taken as their most likely values from Ta-ble 2.To connect with a model for a disk eclipsing the in-ner binary of KH 15D, we approximate the inner and outeredges of the disk as two rings with different orbital elements { r k , I k , Ω k } , with k = L , T for the leading and tailing edges ofthe disk, respectively. Although each ring has a different r k ,we assume the rings precesses about ˆ l with the same globaldisk precession frequency ν d . MNRAS000
Data (blue points) and fitted model (black lines) displayed in Figure 3, folded over the binary orbital period, prior to theyear 2005. The timespan over which the data and model are folded over is displayed in each figure. Our model reproduces the changingmorphology of the light-curve of KH 15D well. orbital elements { a , e , ω, I , Ω } , the position of the ring in theframe of the observer is (cid:169)(cid:173)(cid:171) x y z (cid:170)(cid:174)(cid:172) r , obs = R Z ( Ω ) R X ( I ) R Z ( ω ) R Z ( Ω r ) R X ( I r ) (cid:169)(cid:173)(cid:171) x r y r (cid:170)(cid:174)(cid:172) , (7)where ( x r , y r ) = r ( cos ϕ, sin ϕ ) parameterizes the ( X , Y ) coor-dinates of the ring in the frame where ˆ z = ˆ l , and R X [ β ] ( R Z [ γ ] ) denote rotations along the X ( Z )-axis by angles β ( γ ). As in Winn et al. (2006), we choose the reference planeof the observer so the X -axis points along the binary line-of-nodes (so Ω = ). Also, because our MCMC model highly favors a nearly edge-on orbit (Table 2), we assume I (cid:39) ◦ for simplicity for the rest of this section. All other orbitalparameters are taken as their most likely values from Ta-ble 2.To connect with a model for a disk eclipsing the in-ner binary of KH 15D, we approximate the inner and outeredges of the disk as two rings with different orbital elements { r k , I k , Ω k } , with k = L , T for the leading and tailing edges ofthe disk, respectively. Although each ring has a different r k ,we assume the rings precesses about ˆ l with the same globaldisk precession frequency ν d . MNRAS000 , 1–15 (2020)
Poon, Zanazzi, & Zhu
Figure 5.
Same as Figure 4, except for data and model fits after the year 2005. We also display model predictions from the years 2019to 2021.
To connect the geometry of a disk eclipsing the binary ofKH 15D with the occulting edges of the light-curve model inSection 2, we approximate an eclipsing ring by a line drawntangent to the ring at the location where the ring intersectsthe X -axis of the system (see Fig. 6). The angle between thetangent line and X -axis θ k , as well as the Y -intercept Y k , ofring k are then given by θ k [ I k ( t ) , Ω k ( t )] = tan − (cid:20) tan I k sin ( ω + Ω k ) (cid:21) , (8) Y k [ I k ( t ) , Ω k ( t )] = − r k tan I k tan ( ω + Ω k ) . (9)A successful model of the circumbinary disk of KH 15D would give values for θ k and Y k which match the MCMCfits for θ L , θ T , Y L , and Y T , from Table 2. Before presenting an example warped disk geometry whichmatches the light-curve model fit, we discuss how the warpeddisk geometry can be constrained by the MCMC fits of Sec-tion 2.3. To do this, we simplify equations (8)-(9), and deriveorder-of-magnitude estimates for all disk quantities. Becausethe pericenter direction of the binary is nearly perpendicu-lar to the observer ( ω (cid:28) ), the disk annuli longitude of MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Table 2.
Model fits to photometric and radial velocity datafor the KH 15D system, taking α = . . Orbital parameters { P , e , I , ω, T p } are constrained using photometry and radial ve-locity data, while the other parameters are constrained using pho-tometry alone. Parameter Our Fit P [days] . + . − . e . + . − . I [deg] . + . − . ω [deg] . + . − . T p [JD] - 2,452,350 . + . − . L B / L A . + . − . (cid:15) . + . − . (cid:15) . + . − . ξ . + . − . ξ . + . − . t . + . − . t . + . − . t . + . − . θ L ( t ) [deg] − . + . − . θ T ( t ) [deg] − . + . − . (cid:219) θ L [rad/year] . + . − . (cid:219) θ L [rad/year] . + . − . (cid:219) θ T [rad/year] − . + . − . χ χ t . ± . t . ± . t . ± . t . ± . t . ± . t . ± . t ± Y T ( t ) [au] -0.05678 Y L ( t ) [au] 0.06966 Y T ( t ) [au] -0.02196 (cid:219) Y L ( t ) [au/year] 0.003362 (cid:219) Y T ( t ) [au/year] 0.001681 Predicted by the free parameters. Best-fit value, we do not calculate theerrors implied by t i measurements. ascending nodes satisfy Ω k ≈ π / during transit. The diskinclination is also nearly aligned with the orbital plane of thebinary ( | I k | (cid:28) ). Also, because the binary pericenter direc-tion is nearly perpendicular to the observer, the inclinationnutations should be near a local minimum (e.g. Farago &Laskar 2010; Zanazzi & Lai 2018), so (cid:219) I k ≈ . The nodalregression rate of the rings should be of order (cid:219) Ω k ≈ − ν k ,where ν k is the nearly constant nodal precession rate of ring k . Defining δ Ω k ≡ Ω k − π / , and assuming | ω | , | I k | , and | δ Ω k | (cid:28) , equations (8)-(9) can be shown to reduce to θ k ≈ I k (10) (cid:219) θ k ≈ −( ω + δ Ω k ) I k ν k , (11) Y k ≈ r k I k ( ω + δ Ω k ) , (12) (cid:219) Y k ≈ − r k I k ν k . (13)From this, we see the increase of Y L and Y T is primarily dueto nodal regression from the rings. The evolution of θ L and θ T is primarily due to the curvature of the ring, as it nodally Y(t) Y X θ (t) Star B Star A (x r, , y r ) Test particle ringOcculting Edge
Figure 6.
Our interpretation for the leading/tailing edges ofthe opaque screen in our photometry model. The leading/tailingedges of the screen are from the inner or outer disk truncationradii. The occulting disk edge is approximated by a straight line,drawn tangent to the intersection of the ring with the X -axis ofthe coordinate system. The θ = θ k and Y = Y k values of theleading/tailing edge are defined similarly as in Figure 2. Becausethe binary is nearly edge-on, the straight-line approximation isexcellent. precesses in front of the orbit of the binary (see Silvia & Agol2008 for further discussion). Most interestingly, the MCMCconstraints on θ L and θ T directly translate to constraints onthe ring inclinations I L and I T .Assuming the disk precesses rigidly ( ν L ≈ ν T ), one canthen constrain the disk radial extent. Equation (13) leads to r L r T ≈ θ T θ L (cid:219) Y L (cid:219) Y T = . (cid:18) θ T − ◦ (cid:19) (cid:18) − ◦ θ L (cid:19) (cid:18) . α (cid:19) . (14)Because the values of α which fit the data are of order unity( . (cid:46) α (cid:46) ), we can be confident that the leading edgeof the screen eclipsing the binary of KH 15D is the diskinner truncation radius, while the tailing edge is the outertruncation radius ( r L (cid:46) r T ).Moreover, because (cid:219) θ k , Y k , and (cid:219) Y k are all known, one canget unique solutions for r k , δ Ω k , and ν k . Starting with ν k ,equations (11)-(13) can be re-arranged to give ν k ≈ (cid:18) (cid:219) θ k (cid:219) Y k Y k (cid:19) / . (15)Evaluating estimate (15) at t = t , we find ν L ∼ .
026 yr − and ν T ∼ . − , which are consistent with one anotherwithin a factor of a few. Similarly, equation (13) can besolved for r k : r k ≈ − θ k (cid:18) Y k (cid:219) Y k (cid:219) θ k (cid:19) / , (16)which gives r L ∼ .
54 au and r T ∼ . at t = t for ourmodel. Last, either equation (11) or (12) can be solved for δ Ω k .Although these disk parameter estimates are far from MNRAS000
54 au and r T ∼ . at t = t for ourmodel. Last, either equation (11) or (12) can be solved for δ Ω k .Although these disk parameter estimates are far from MNRAS000 , 1–15 (2020) Poon, Zanazzi, & Zhu
Parameter Example Value ν d [ yr − ] 0.01 r L [AU] 0.5 r T [AU] 2.0 I L ( t ) [deg] -13 I T ( t ) [deg] -6 Ω L ( t ) [deg] 100 Ω T ( t ) [deg] 115 Table 3.
Parameter values for our dynamical model of a warpeddisk precessing around an eccentric binary in the KH 15D system.Disk inclinations are relative to the binary orbital plane, and disklongitude of ascending nodes are relative to the binary pericenterdirection. See text for definitions and discussion. unique, they provide constraints on the properties of the cir-cumbinary disk within the KH 15D system. We can stronglyconclude the disk-binary mutual inclination I KH 15D in theKH 15D system lies in the range ◦ (cid:46) I KH 15D (cid:46) ◦ , withthe disk inner edge more highly inclined than the outer edge(because I L (cid:38) I T ). The leading edge of the opaque screencrossing the binary orbit is from the disk inner edge, whichis located at a radius r L (cid:46) , while the tailing outer diskedge is located at r T ∼ few au . As we saw in the previous section, for a unique match tothe phenomenological parameters { θ k , (cid:219) θ k , Y k , (cid:219) Y k } to a pre-cessing, inclined ring annulus, we require the ring param-eters { r k , I k , Ω k , ν k } . However, for a protoplanetary disk toexist over many dynamical times, it must precess rigidly( ν L = ν T ), decreasing the number of free parameters in onering. Therefore, our dynamical model is over-determined byour phenomenological model. To get accurate constraints onthe warped disk itself using photometry, a light-curve modelmust be developed whose free parameters are directly relatedto the warped disk properties (disk inclination, warp, twist,precession frequency, etc.), rather than indirectly througha phenomenological model. The goal of this section is notto provide stringent constraints on the disk itself, but topresent an example warped disk which gives gross light-curvefeatures consistent with the MCMC light-curve fits.Motivated by the estimated leading and tailing esti-mates in the previous subsection, we experiment with thewarped disk orbital parameters and global precession fre-quency, to find a disk whose properties match the MCMC fit-ted parameters. Table 3 presents example model parametersfor a dynamically evolving, warped disk whose features arecompatible with the light-curve fits, with Figure 7 display-ing θ k ( t ) and Y k ( t ) for both (dynamical and MCMC) modelsover the duration of time the leading and tailing eclipseshave been observed. The dynamical model and MCMC fitsmatch one another within a factor of a few, heavily reinforc-ing the idea that the light-curve of KH 15D is caused bya warped, relatively narrow, precessing disk, eclipsing thestarlight of the eccentric binary. In particular, we see thebehavior of the dynamical model matches the θ L ( t ) light-curve fit, reproducing the decrease in (cid:219) θ L before and afterthe year t = t (cid:39) . Figure 7.
Comparing our dynamical warped disk model (Ta-ble 3) with the MCMC fits from our phenomenological photome-try model (Table 2). Although agreement between the two modelscan be improved, the warped disk reproduces the main featuresof the photometry model.
Complimentary constraints on the disk of KH 15D havecome recently from the double-peaked line profile of neutraloxygen emission, assuming the [ OI ] λ emission origi-nates from the surface of the gaseous circumbinary disk ofKH 15D (Fang et al. 2019). This line profile was used toconstrain the disk radial extent, as well as the disk surfacedensity profile. Fang et al. (2019) found an inner disk radiusof r in ≈ .
57 au , an outer radius r out ≈ . , and surface den-sity profile Σ ∝ r − . . The inner and outer radii are roughlyconsistent with our r L = r in and r T = r out values constrainedby our crude dynamical fit to the photometry of KH 15D(Table 3). We note that the outer edge of a protoplanetarydisk gas and dust radius may differ, due to radial drift ofthe dust (e.g. Weidenschilling 1977; Takeuchi & Lin 2002;Birnstiel & Andrews 2014; Powell et al. 2017; Rosotti et al.2019). Indeed, molecular line and continuum emission havebeen shown to extend to different radii around young stel-lar objects (e.g. Pani´c et al. 2009; Andrews et al. 2012; deGregorio-Monsalvo et al. 2013; Ansdell et al. 2018; Facchini MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Figure 8.
Global disk precession frequnecy ν d (eq. 17) as a func-tion of the surface density powerlaw index p ( Σ ∝ r p ), for theKH 15D system parameters, assuming r in = . with the r out values indicated. Dotted green line shows our dynamical modelvalue of ν d = .
01 yr − . Depending on the disk radial extent, ameasurement of ν d translates to a constraint on p . et al. 2019), showing gas and dust in protoplanetary disksoften extend out to different radii (sometimes differing by asmuch as a factor of ∼ ). The dynamics of dust in a precess-ing circumbinary disk can also be non-trivial (Poblete et al.2019; Aly & Lodato 2020).Our warped disk model can also constrain the disk sur-face density profile Σ ∝ r p , because the distribution of masswithin the disk affects the torque exerted on the disk by thebinary, modifying the disk precession frequency ν d . Assum-ing a nearly-flat disk which is driven into rigid-body preces-sion about the binary, ν d can be shown to be (e.g. Lodato& Facchini 2013; Foucart & Lai 2014; Zanazzi & Lai 2018;Lubow & Martin 2018) ν d = (cid:18) / + p − p (cid:19) (cid:20) − ( r out / r in ) p − ( r out / r in ) / + p − (cid:21) µ M t (cid:18) ar in (cid:19) (cid:115) GM t r . (17)Figure 8 plots the ν d value given by equation (17) as a func-tion of p . Depending on the disk outer radius, we clearly seea measurement of ν d can constrain the p value of the disk.The model parameters from Table 3 support p ∼ , whichdiffers substantially from the Fang et al. (2019) constraintof p ≈ − . . Further photometric modelling is required tosee if the p value implied by the disk precession frequencydiffers from that constrained by the disk OI emission. In §
3, we showed how the photometry of KH 15D could beexplained by a precessing circumbinary disk, in agreementwith the results of other works (Winn et al. 2004, 2006;Chiang & Murray-Clay 2004; Silvia & Agol 2008). The pa-rameters constrained by the photometry of the system arelisted in Table 3. Although the fit of the dynamical model to the photometry is crude, we argue the basic conclusions onthe parameters of the system are unlikely to differ by morethan a factor of a few, and comprises some of the first con-straints on small warps within protoplanetary disks. Thissection connects the constraints of our dynamical model totheories describing warp propagation in accretion disks, aswell as speculation on the long-term evolution of the system.We also discuss predictions from our model, as well as futuremodelling efforts.
Our dynamical model requires a non-zero warp ( ∆ I = I T − I L )and twist ( ∆Ω = Ω T − Ω L ) to cause the complex series ofeclipses seen in the KH 15D system. These warps and twistsarise from the disk resisting the differential nodal precessioninduced by the specific torque of the binary | T bin | ∼ r n ν | ¯ I | , (18)where n = (cid:112) GM t / r is the rings orbital frequency, ν is thecharacteristic nodal precession frequency induced on the diskfrom the binary (eq. 6), and ¯ I is the characteristic “average”inclination of the disk. One way to balance the torque fromthe binary is by thermal pressure between ringlets, whichhas an internal torque of order (e.g. Ogilvie 1999; Chiang &Culter 2003; Chiang & Murray-Clay 2004) | T press | ∼ c | ∆ I | , (19)where c s = hrn is the ring sound-speed, while h is the aspectratio of the disk. Assuming torque balance ( | T bin | ≈ | T press | )allows us to estimate the warp which may develop under theresisting influence of thermal pressure: (cid:12)(cid:12)(cid:12)(cid:12) ∆ I ¯ I (cid:12)(cid:12)(cid:12)(cid:12) press ∼ r n ν c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = ¯ r = (cid:18) . h (cid:19) (cid:18) µ .
33 M (cid:12) (cid:19) (cid:18) .
32 M (cid:12) M t (cid:19) (cid:16) a .
29 au (cid:17) (cid:18) . r (cid:19) , (20)where ¯ r is some characteristic radius within the disk. Clearly,this warp is quite large.However, in nearly-inviscid (Shakura-Sunyaev parame-ter α (cid:46) h ) disks with the radial-epicyclic frequency satisfy-ing κ ≈ n , the near-resonant propagation of bending wavesacross the disk can amplify the strength of the hydrody-namical torque by a factor (Papaloizou & Lin 1995; Lubow& Ogilvie 2000; Ogilvie 2006) | T bw | ∼ | ˜ κ | | T press | , (21)where ˜ κ ≡ ( κ − n )/( n ) is a dimensionless quantity relatedto the apsidal precession rate. Because the secular apsidalprecession rate is | ˜ κ | ∼ ν / n for circumbinary disks (Miranda& Lai 2015), torque balance ( | T bw | ≈ | T bin | ) gives (cid:12)(cid:12)(cid:12)(cid:12) ∆ I ¯ I (cid:12)(cid:12)(cid:12)(cid:12) bw ∼ r ν c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = ¯ r = . (cid:18) . h (cid:19) (cid:18) µ .
33 M (cid:12) (cid:19) (cid:18) .
32 M (cid:12) M t (cid:19) (cid:16) a .
29 au (cid:17) (cid:18) . r (cid:19) . (22) MNRAS000
29 au (cid:17) (cid:18) . r (cid:19) . (22) MNRAS000 , 1–15 (2020) Poon, Zanazzi, & Zhu
This estimate is much closer to the | ∆ I |/ ¯ I ∼ values impliedby our dynamical model, and lies in agreement with moredetailed calculations of warp propagation in protoplanetarydisks (Lodato & Facchini 2013; Foucart & Lai 2013, 2014;Zanazzi & Lai 2018; Lubow & Martin 2018).Disk self-gravity can also resist differential nodal pre-cession from the binary. Mutually-misaligned ringlets expe-rience specific mutual internal torques of order (Chiang &Culter 2003; Chiang & Murray-Clay 2004; Tremaine & Davis2014; Zanazzi & Lai 2017; Batygin 2018) | T sg | ∼ G Σ rh | ∆ I | ∼ GM d hr | ∆ I | , (23)assuming the disk mass M d ∼ r Σ , and the additional factorof h − arises from the enhancement of the mutual gravita-tional attraction between ringlets when the disk is verticallythin (Batygin 2018). Torque balance ( | T sg | ≈ | T bin | ) leads towarps of order (cid:12)(cid:12)(cid:12)(cid:12) ∆ I ¯ I (cid:12)(cid:12)(cid:12)(cid:12) sg ∼ r n ν hGM d (cid:12)(cid:12)(cid:12)(cid:12) r = ¯ r = . (cid:18) h . (cid:19) (cid:18) µ .
33 M (cid:12) (cid:19) (cid:18) . Jup M d (cid:19) (cid:16) a .
29 au (cid:17) (cid:18) . r (cid:19) . (24)Even after assuming the upper limit on the total (gas anddust) disk mass inferred by ALMA observations (Aronowet al. 2018), self-gravity is typically not as effective as bend-ing waves at enforcing coplanarity between ringlets. How-ever, a massive disk ( M d ∼ Jup ) can give warps compera-ble to those inferred by our dynamical KH 15D disk model.The direction of the warp ( ∆ I positive or negative)has also been argued to encode information on the internalforces/torques enforcing disk coplanarity. Chiang & Murray-Clay (2004) argued thermal pressure predicts ∆ I < , whileself-gravity predicts ∆ I > . More detailed calculations sup-port the prediction a disk should relax to a ∆ I > pro-file under the influence of disk self-gravity (Batygin 2012,2018; Zanazzi & Lai 2017). But calculations taking intoaccount the resonant propagation of bending waves also predict ∆ I > (e.g. Facchini et al. 2013; Foucart & Lai2014; Zanazzi & Lai 2018; Lubow & Martin 2018). Hydro-dynamical simulations of protoplanetary disks (neglectingself-gravity) find conflicting results, with ∆ I > and ∆ I < at different times, primarily because the simulations usu-ally cannot be run long enough for the system to relax toa smoothly-evolving warp profile (e.g. Facchini et al. 2013;Martin & Lubow 2017, 2018; Smallwood et al. 2019, 2020;Moody et al. 2019). We note the disk may never relax toa steady-state. Simulations which accurately calculate howthe binary interacts with a tidally-truncated circumbinarydisk find highly-dynamical inner disk edges for disks orbit-ing eccentric binaries (e.g. Miranda et al. 2017; Mu˜noz et al.2019, 2020; Franchini et al. 2019). Because resonant Lind-blad torques often truncate disks (e.g. Artymowicz & Lubow1994; Lubow et al. 2015; Miranda & Lai 2015), which mayalso excite disk tilts (Borderies et al. 1984; Lubow 1992;Zhang & Lai 2006), it is not unreasonable to say a real cir-cumbinary disk may never relax to a steady-state inclinationprofile. We conclude that bending-wave propagation is themain internal force enforcing rigid precession of the disk ofKH 15D, despite the conflicting predictions for the sign of ∆ I . A small viscosity in a circumbinary disk also leads to anon-zero twist, due to the azimuthal shear induced by differ-ential nodal precession. The magnitude of the torque resist-ing nodal shear is (Papaloizou & Pringle 1983; Papaloizou& Lin 1995; Ogilvie 1999; Lubow & Ogilvie 2000) | T visc | ∼ α c | ∆Ω | , (25)assuming an isotropic kinematic viscosity ( ν = α p /[ ρ n ] ).The α − (rather than α + ) dependence in equation (25) isfrom near-resonant forcing of radial and azimuthal pertur-bations (since κ ≈ n ), which are damped only by viscosity(Papaloizou & Lin 1995; Lubow & Ogilvie 2000; Lodato &Pringle 2007). Viscosity leads to twists of order (assuming | T visc | ≈ | T bin | ) (cid:12)(cid:12)(cid:12)(cid:12) ∆Ω ¯ I (cid:12)(cid:12)(cid:12)(cid:12) visc ∼ α r n ν c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = ¯ r = . (cid:16) α . (cid:17) (cid:18) . h (cid:19) (cid:18) µ .
33 M (cid:12) (cid:19) × (cid:18) .
32 M (cid:12) M t (cid:19) (cid:16) a .
29 au (cid:17) (cid:18) . r (cid:19) . (26)More detailed calculations typically give positive ∆Ω values abit larger in circumbinary disks (Foucart & Lai 2014; Zanazzi& Lai 2018), in agreement with our dynamical model. Al-though observations frequently infer α values much lowerthan − (e.g. Hughes et al. 2011; Flaherty et al. 2015;Teague et al. 2016; Rafikov 2017; Ansdell et al. 2018), thelarge warp in this disk can excite parametric instabilities, en-hancing the viscous dissipation rate in the disk (Goodman1993; Ryu & Goodman 1994; Gammie et al. 2000; Ogilvie& Latter 2013; Paardekooper & Ogilvie 2019).We conclude the disk warp and twist implied by ourmodel lie in accord with hydrodynamical theories of warpedaccretion disks. Recently, Czekala et al. (2019) showed circumbinary disks(both gas and debris) have higher inclinations when orbitingeccentric binaries. Figure 9 displays the disk inclinations an-alyzed in Czekala et al. (2019), alongside our constraints forthe inclination of KH 15D, which we take directly from ourphotometric fits ( | θ T [ t ]| (cid:46) | I KH 15D | (cid:46) | θ L [ t ]| , see § I crit = cos − (cid:115) e + e , (27)which is a necessary (but not sufficient) condition for thedisk-binary inclination to evolve to ◦ (polar alignment).From Figure 9, because | I crit | > | I KH 15D | , we can be confi-dent the disk will not polar align, and will eventually alignwith the orbital plane of the binary (without any other mech-anisms exciting the disk inclination).We can also estimate the timescale over which the diskinclination evolves. Because a non-zero twist ∆Ω exerts abackreaction torque on the disk from the binary, the disk isdriven into alignment (or polar alignment) over the timescale MNRAS , 1–15 (2020) onstraining the tilt in KH 15D M u t u a l I n c li n a t i o n [ d e g ] KH 15DProtoplanetary disksDebris disks
Figure 9.
The mutual inclination between the disk and binaryorbital planes in the KH 15D system (red), plotted alongside cir-cumbinary disk inclinations for protoplanetary (orange) and de-bris (green) disks (Czekala et al. 2019), as a function of binaryeccentricity. The dashed blue lines plots the critical inclination(eq. 27). The black dotted lines connect degenerate solutions forHD 142527, SR 24N and GG Tau Aa-Ab. The triangle representsthe lower limit for R CrA. The disk in KH 15D will align (notpolar align) with the orbital plane of the binary. (Foucart & Lai 2014; Zanazzi & Lai 2018) γ evol ∼ ν d | ∆Ω | . (28)Inserting equation (26) into equation (28) gives the often-quoted “Bate timescale” (Bate et al. 2000). However, because ν d and ∆Ω are both determined by our dynamical model,we can actually estimate γ evol using observationally inferredparameters : γ evol ∼ . × − (cid:18) ν d .
01 yr − (cid:19) (cid:18) ∆Ω ◦ (cid:19) yr − . (29)Equation (29) implies the disk should align with the orbitalplane of the binary in less than ∼ years. Although secularinteractions can keep the disk misaligned with the eccentricorbital plane of the binary over timescales a few times longerthan estimate (29) (Zanazzi & Lai 2018; Smallwood et al.2019), this is much shorter than the ∼ year lifetimes oftypical protoplanetary disks (e.g. Haisch et al. 2001). Eitherwe are observing KH 15D while it is still very young, oradditional mechanisms are exciting the disk inclination. The most immediate consequences are predictions for futurelight-curve behavior from our photometry model fits (Figs. 3& 5, Table 2). Current I -band measurements should showthe light from star A slowly being revealed by the tailingedge (since t ≈ ). By the year ∼ ∼ (cid:219) θ T , our dynamical model predicts that the fitcan be further improved if the change of (cid:219) θ T with time isincorporated (eq. 11).We are able to make an explicit connection between the phenominological model ( § § { r k , I k , Ω k } and global disk preces-sion frequency ν d , rather than parameters describing the lo-cations and orientations of the leading and trailing edges { θ k , (cid:219) θ k , Y k , (cid:219) Y k } , to fit the light-curve of KH 15D. This exer-cise should yield parameters consistent with those listed inTable 3 within a factor of a few.Our folded light-curves (Figs 4-5) show the leading edgeis well fit by a sharp edge, whereas the poor fit for the tailingedge imply it is clumpy/puffy, in agreement with the find-ings of Garc´ıa Soto et al. (2020). The sharp inner edge islikely due to tidal truncation by the torque from the binary.Calculations and hydrodynamical simulations suggest thatthe radius at which the binary truncates the disk is ∼ a ≈ . au) (e.g. Miranda & Lai2015), lying close to the inner radius value of our dynamicalmodel (Table 3). However, it remains unclear why the diskis so compact ( r out (cid:46) ), and the possibility still exists thedisk outer edge is truncated by a planet (Chiang & Murray-Clay 2004). Future modelling of how dust scattering and thefinite optical depth of the disk can create a “fuzzier” outeredge would be of interest (Chiang & Murray-Clay 2004; Sil-via & Agol 2008). In this work, we have developed a circumbinary disk modelthat explains the photometric variability of KH 15D span-ning more than 60 years. From this model, we are able toconstrain the disk annular extent, inclination, orientationwith respect to the binary pericenter direction, warp pro-file, precession frequency, and even surface density profile.The fits of our phenomenological model to fit the photom-etry of KH 15D are displayed in Table 2, with parametersof a warped disk which are consistent with the phenomeno-logical model constraints listed in Table 3. Although strictconstraints on the warped disk remain elusive, we can beconfident about the following features of the disk: • The beginning of the dips/eclipses in KH 15D are dueto the disk inner edge slowly covering the binary, while thecurrently observed slow reversal of the dipping behavior inKH 15D is due to the disk outer edge slowly revealing thebinary. The inner edge has a radius r in (cid:46) , while theouter edge has a radius r out ∼ few au . • The disk inner edge is more inclined to the orbital planeof the binary than the disk outer edge. Both inner and outerdisk inclinations are less than ∼ ◦ , but greater than ∼ ◦ ,with a difference of order ∼ ◦ . • The disk inner and outer longitude of ascending nodesdiffer by ∼ ◦ .These constraints are consistent with hydrodynamical theo-ries of warped accretion disks, resisting differential nodalprecession from the gravitational torque from the binary( § P prec ∼ π / ν d ∼
600 years , but this constraint is sensitiveto the model fit of KH 15D. We can be very confident, how-ever, that the timescale over which the disk of KH 15D aligns
MNRAS000
MNRAS000 , 1–15 (2020) Poon, Zanazzi, & Zhu with the orbital plane of the binary is much shorter thanthe lifetime of the disk (eq. 29), suggesting that additionalmechanisms are exciting the disk tilt.
ACKNOWLEDGEMENTS
We thank Eugene Chiang, Ruth Murray-Clay, and NormMurray for useful discussions, and Jessica Speedie for help indeveloping the circumbinary disk visualization in Figure 1.MP has been supported by undergraduate research fellow-ships from the Centre for Planetary Sciences and CanadianInstitute for Theoretical Astrophysics (CITA). JZ is sup-ported by a CITA postdoctoral fellowship. WZ is supportedby the Beatrice and Vincent Tremaine Fellowship at CITA.
REFERENCES
Aly H., Lodato G., 2020, MNRAS, 492, 3306Aly H., Dehnen W., Nixon C., King A., 2015, MNRAS, 449, 65Andrews S. M., Czekala I., Wilner D. J., Espaillat C., DullemondC. P., Hughes A. M., 2010, ApJ, 710, 462Andrews S. M., et al., 2012, ApJ, 744, 162Ansdell M., et al., 2018, ApJ, 859, 21Armstrong D. J., Osborn H. P., Brown D. J. A., Faedi F., G´omezMaqueo Chew Y., Martin D. V., Pollacco D., Udry S., 2014,MNRAS, 444, 1873Aronow R. A., Herbst W., Hughes A. M., Wilner D. J., WinnJ. N., 2018, AJ, 155, 47Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651Bate M. R., 2012, MNRAS, 419, 3115Bate M. R., 2018, MNRAS, 475, 5618Bate M. R., Bonnell I. A., Clarke C. J., Lubow S. H., OgilvieG. I., Pringle J. E., Tout C. A., 2000, MNRAS, 317, 773Bate M. R., Bonnell I. A., Bromm V., 2002, MNRAS, 336, 705Batygin K., 2012, Nature, 491, 418Batygin K., 2018, MNRAS, 475, 5070Bi J., et al., 2020, ApJ, 895, L18Birnstiel T., Andrews S. M., 2014, ApJ, 780, 153Borderies N., Goldreich P., Tremaine S., 1984, ApJ, 284, 429Brinch C., Jørgensen J. K., Hogerheijde M. R., Nelson R. P.,Gressel O., 2016, ApJ, 830, L16Capelo H. L., Herbst W., Leggett S. K., Hamilton C. M., JohnsonJ. A., 2012, ApJ, 757, L18Chen C., Franchini A., Lubow S. H., Martin R. G., 2019, MNRAS,490, 5634Chen C., Lubow S. H., Martin R. G., 2020, MNRAS, 494, 4645Chiang E. I., Culter C. J., 2003, ApJ, 599, 675Chiang E. I., Murray-Clay R. A., 2004, ApJ, 607, 913Czekala I., Andrews S. M., Jensen E. L. N., Stassun K. G., TorresG., Wilner D. J., 2015, ApJ, 806, 154Czekala I., Andrews S. M., Torres G., Jensen E. L. N., StassunK. G., Wilner D. J., Latham D. W., 2016, ApJ, 818, 156Czekala I., et al., 2017, ApJ, 851, 132Czekala I., Chiang E., Andrews S. M., Jensen E. L. N., Torres G.,Wilner D. J., Stassun K. G., Macintosh B., 2019, ApJ, 883,22Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656Doyle L. R., Deeg H. J., 2018, The Way to Circumbinary Planets.p. 115, doi:10.1007/978-3-319-55333-7 115Facchini S., Lodato G., Price D. J., 2013, MNRAS, 433, 2142Facchini S., Juh´asz A., Lodato G., 2018, MNRAS, 473, 4459Facchini S., et al., 2019, A&A, 626, L2Fang M., Pascucci I., Kim J. S., Edwards S., 2019, ApJ, 879, L10Farago F., Laskar J., 2010, MNRAS, 401, 1189 Flaherty K. M., Hughes A. M., Rosenfeld K. A., Andrews S. M.,Chiang E., Simon J. B., Kerzner S., Wilner D. J., 2015, ApJ,813, 99Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Foucart F., Lai D., 2013, ApJ, 764, 106Foucart F., Lai D., 2014, MNRAS, 445, 1731Franchini A., Lubow S. H., Martin R. G., 2019, ApJ, 880, L18Fulton B. J., Petigura E. A., Blunt S., Sinukoff E., 2018, PASP,130, 044504Gammie C. F., Goodman J., Ogilvie G. I., 2000, MNRAS, 318,1005Garc´ıa Soto A., Ali A., Newmark A., Herbst W., Windemuth D.,Winn J. N., 2020, AJ, 159, 135Giuppone C. A., Cuello N., 2019, in Journal of Physics Confer-ence Series. p. 012023 ( arXiv:1907.08180 ), doi:10.1088/1742-6596/1365/1/012023Goodman J., 1993, ApJ, 406, 596Haisch Karl E. J., Lada E. A., Lada C. J., 2001, ApJ, 553, L153Hamilton C. M., Herbst W., Mundt R., Bailer-Jones C. A. L.,Johns-Krull C. M., 2003, ApJ, 591, L45Hamilton C. M., et al., 2005, AJ, 130, 1896Hughes A. M., Wilner D. J., Andrews S. M., Qi C., HogerheijdeM. R., 2011, ApJ, 727, 85Johnson J. A., Winn J. N., 2004, AJ, 127, 2344Johnson J. A., Marcy G. W., Hamilton C. M., Herbst W., Johns-Krull C. M., 2004, AJ, 128, 1265Kearns K. E., Herbst W., 1998, ApJ, 116, 261Kennedy G. M., et al., 2012a, MNRAS, 421, 2264Kennedy G. M., Wyatt M. C., Sibthorpe B., Phillips N. M.,Matthews B. C., Greaves J. S., 2012b, MNRAS, 426, 2115Kennedy G. M., et al., 2019, Nature Astronomy, 3, 230Kraus S., et al., 2020, arXiv e-prints, p. arXiv:2004.01204Larwood J. D., Papaloizou J. C. B., 1997, MNRAS, 285, 288Li G., Holman M. J., Tao M., 2016, ApJ, 831, 96Lodato G., Facchini S., 2013, MNRAS, 433, 2157Lodato G., Pringle J. E., 2007, MNRAS, 381, 1287Lubow S. H., 1992, ApJ, 398, 525Lubow S. H., Martin R. G., 2018, MNRAS, 473, 3733Lubow S. H., Ogilvie G. I., 2000, ApJ, 538, 326Lubow S. H., Martin R. G., Nixon C., 2015, ApJ, 800, 96Maffei P., Ciprini S., Tosti G., 2005, MNRAS, 357, 1059Marino S., Perez S., Casassus S., 2015, ApJ, 798, L44Martin R. G., Lubow S. H., 2017, ApJ, 835, L28Martin R. G., Lubow S. H., 2018, MNRAS, 479, 1297Martin D. V., Triaud A. H. M. J., 2014, A&A, 570, A91McLaughlin D. B., 1924, ApJ, 60, 22Miranda R., Lai D., 2015, MNRAS, 452, 2396Miranda R., Mu˜noz D. J., Lai D., 2017, MNRAS, 466, 1170Moody M. S. L., Shi J.-M., Stone J. M., 2019, ApJ, 875, 66Mu˜noz D. J., Miranda R., Lai D., 2019, ApJ, 871, 84Mu˜noz D. J., Lai D., Kratter K., Mirand a R., 2020, ApJ, 889,114Nixon C., King A., Price D., 2013, MNRAS, 434, 1946Ogilvie G. I., 1999, MNRAS, 304, 557Ogilvie G. I., 2006, MNRAS, 365, 977Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2420Paardekooper S.-J., Ogilvie G. I., 2019, MNRAS, 483, 3738Pani´c O., Hogerheijde M. R., Wilner D., Qi C., 2009, A&A, 501,269Papaloizou J. C. B., Lin D. N. C., 1995, ApJ, 438, 841Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181Poblete P. P., Cuello N., Cuadra J., 2019, MNRAS, 489, 2204Poblete P. P., Calcino J., Cuello N., Mac´ıas E., Ribas ´A., PriceD. J., Cuadra J., Pinte C., 2020, MNRAS, 496, 2362Powell D., Murray-Clay R., Schlichting H. E., 2017, ApJ, 840, 93Price D. J., et al., 2018, MNRAS, 477, 1270Rafikov R. R., 2017, ApJ, 837, 163 MNRAS , 1–15 (2020) onstraining the tilt in KH 15D Rosenfeld K. A., Andrews S. M., Wilner D. J., Stempels H. C.,2012, ApJ, 759, 119Rosotti G. P., Tazzari M., Booth R. A., Testi L., Lodato G.,Clarke C., 2019, MNRAS, 486, 4829Rossiter R. A., 1924, ApJ, 60, 15Ru´ız-Rodr´ıguez D., Kastner J. H., Dong R., Principe D. A., An-drews S. M., Wilner D. J., 2019, AJ, 157, 237Ryu D., Goodman J., 1994, ApJ, 422, 269Silvia D. W., Agol E., 2008, ApJ, 681, 1377Smallwood J. L., Lubow S. H., Franchini A., Martin R. G., 2019,MNRAS, 486, 2919Smallwood J. L., Franchini A., Chen C., Becerril E., Lubow S. H.,Yang C.-C., Martin R. G., 2020, MNRAS, 494, 487Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344Teague R., et al., 2016, A&A, 592, A49Tremaine S., Davis S. W., 2014, MNRAS, 441, 1408Weidenschilling S. J., 1977, MNRAS, 180, 57Welsh W. F., Orosz J. A., 2018, Two Suns in the Sky: TheKepler Circumbinary Planets. p. 34, doi:10.1007/978-3-319-55333-7 34Winn J. N., Holman M. J., Johnson J. A., Stanek K. Z., GarnavichP. M., 2004, ApJ, 603, L45Winn J. N., Hamilton C. M., Herbst W. J., Hoffman J. L., HolmanM. J., Johnson J. A., Kuchner M. J., 2006, ApJ, 644, 510Zanazzi J. J., Lai D., 2017, MNRAS, 464, 3945Zanazzi J. J., Lai D., 2018, MNRAS, 473, 603Zhang Z., Fabrycky D. C., 2019, ApJ, 879, 92Zhang H., Lai D., 2006, MNRAS, 368, 917de Gregorio-Monsalvo I., et al., 2013, A&A, 557, A133
APPENDIX A: MODEL PARAMETERS
We display all MCMC parameter fits for our new photomet-ric model, for various α , as described in Section 2.2. Recall α is the ratio of the tail edge velocity over the lead edge veloc-ity along the vertical axis of our line of sight. Small α testsfor narrow disk models whereas large α tests for extendeddisk models. We do not report model parameters for α = . since the MCMC does not converge. Upper and lower errorbars indicate a σ confidence interval. Model parameters aredescribed in Table 1. APPENDIX B: MCMC CORNER PLOTS
We display the corner plots to our best fit model ( α = . )with best fit values listed in Table 2. We remove the first17,500 of 20,000 total steps as burn-in, and plot the poste-rior distribution. The apparent degeneracy with ξ and L B appears in many of the MCMC fits. This is likely due tosome subtleties of the halo model, yet they do not affectthe quality of the photometric fits. Because we are primar-ily interested in constraints on the properties describing theascent of the leading and tailing screens { θ k , (cid:219) θ k , Y k , (cid:219) Y k } , to beconsistent with Winn et al. (2006), we do not modify thehalo model. This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS , 1–15 (2020) Poon, Zanazzi, & Zhu
Table A1.
Same as Table 2, except we vary the value of α .Free parameter α = . α = . α = . α = . P [days] . + . − . . + . − . . + . − . . + . − . e . + . − . . + . − . . + . − . . + . − . I [deg] . + . − . . + . − . . + . − . . + . − . ω [deg] . + . − . . + . − . . + . − . . + . − . T p [JD] - 2,452,350 . + . − . . + . − . . + . − . . + . − . L B / L A . + . − . . + . − . . + . − . . + . − . (cid:15) . + . − . . + . − . . + . − . . + . − . (cid:15) . + . − . . + . − . . + . − . . + . − . ξ . + . − . . + . − . . + . − . . + . − . ξ . + . − . . + . − . . + . − . . + . − . t . + . − . . + . − . . + . − . . + . − . t . + . − . . + . − . . + . − . . + . − . t . + . − . . + . − . . + . − . . + . − . θ L ( t ) [deg] − . + . − . − . + . − . − . + . − . − . + . − . θ T ( t ) [deg] . + . − . − . + . − . − . + . − . − . + . − . (cid:219) θ L [rad/year] . + . − . . + . − . . + . − . . + . − . (cid:219) θ L [rad/year] . + . − . . + . − . . + . − . . + . − . (cid:219) θ T [rad/year] − . + . − . − . + . − . . + . − . . + . − . Fit photometry? No Yes No No χ χ
15 15 18 17 MNRAS , 1–15 (2020) onstraining the tilt in KH 15D T p eωIL B ε ε ξ ξ P e ω I L B ε ε ξ ξ t t θ L1 T p t θ L2 θ T t t t θ L (t ) θ T (t ) θ L1 θ L2 θ T θ L (t ) θ T (t ) P Figure B1.
Two dimensional projection of the posterior probability distribution sampled using MCMC for α = . . Blue solid linesindicate best fit values reported in Table 2 whereas black dashed lines indicate a σ confidence interval.MNRAS000