Table Of Content

MNRAS000,1–18(2016) Preprint11January2017 CompiledusingMNRASLATEXstylefilev3.0 Observational signatures of linear warps in circumbinary discs Attila Juha´sz1(cid:63), Stefano Facchini2† 1Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 2Max-Planck-Institut fu¨r Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany 7 AcceptedXXX.ReceivedYYY;inoriginalformZZZ 1 0 2 ABSTRACT n In recent years an increasing number of observational studies have hinted at the pres- a ence of warps in protoplanetary discs, however a general comprehensive description J of observational diagnostics of warped discs was missing. We performed a series of 0 3D SPH hydrodynamic simulations and combined them with 3D radiative transfer 1 calculations to study the observability of warps in circumbinary discs, whose plane is misaligned with respect to the orbital plane of the central binary. Our numerical ] P hydrodynamic simulations confirm previous analytical results on the dependence of E the warp structure on the viscosity and the initial misalignment between the binary . andthedisc.Tostudytheobservationalsignaturesofwarpswecalculateimagesinthe h continuumatnear-infraredandsub-millimetrewavelengthsandinthepurerotational p transitionofCOinthesub-millimetre.Warpedcircumbinarydiscsshowsurfacebright- - o nessasymmetryinnear-infraredscatteredlightimagesaswellasinopticallythickgas r lines at sub-millimetre wavelengths. The asymmetry is caused by self-shadowing of t s the disc by the inner warped regions, thus the strength of the asymmetry depends on a the strength of the warp. The projected velocity field, derived from line observations, [ showscharacteristicdeviations,twistsandachangeintheslopeoftherotationcurve, 1 from that of an unperturbed disc. In extreme cases even the direction of rotation v appears to change in the disc inwards of a characteristic radius. The strength of the 1 kinematical signatures of warps decreases with increasing inclination. The strength of 1 all warp signatures decreases with decreasing viscosity. 6 2 Key words: accretion, accretion discs — circumstellar matter — protoplanetary 0 discs — hydrodynamics . 1 0 7 1 1 INTRODUCTION Kotze & Charles 2012, and references therein), and some : v cataclysmic variables (summarised by Olech et al. 2009). Protoplanetarydiscsareusuallythoughtofasaxisymmetric i However, very recently, observational features in a few pro- X structures orbiting a single star. However, it is well known toplanetarydiscshavebeeninterpretedassignaturesofdisc that the geometry of such discs can be more complicated. r warping. In particular, the advent of interferometric obser- a In particular, discs in binary systems can be warped. The vations with the Acatama Large Millimeter Array (ALMA) orbitalplaneofthebinaryandtheplaneoftheouterdiscare hasallowedustostartmeasuringprecisegaskinematicsvia likely to be misaligned due to interaction with other stars molecular line profiles, which is one of the most powerful or to accretion of randomly oriented gas (e.g. Bonnell et al. tools to detect disc warping. For example, Rosenfeld et al. 1992;Bateetal.2010).Suchmisalignmentinducesatorque (2012)haveobservedanexcessintheemissionathighveloc- in the disc, leading to a warped 3D structure. ity in the inner disc of TW Hya, which they interpreted as Untilafewyearsago,welackedobservationalevidence a possible inner warp. Another case is HD 100546, where a of warped protoplanetary discs. In binary systems, tilted warphasbeeninvokedtointerpretthetwistedfirstmoment discs had been inferred only around X-ray binaries show- map of CO (3-2) (Pineda et al. 2014), and even earlier as a ing a long term modulation of the periodic light curve (see possibleexplanationofthespiralstructureobservedinscat- teredlight(Quillen2006)andoftheasymmetriclineprofiles fromsingle-dishobservations(Pani´cetal.2010).Averysig- (cid:63) [email protected] nificant case is HD 142527. Marino et al. (2015) have mod- † [email protected] (cid:13)c 2016TheAuthors 2 Juh´asz, Facchini elledthestronglyasymmetricscatteredlightimageofsucha ulations have shown that the warp evolution equations (in protoplanetarydiscwithahighlymisalignedinnerdisc,cast- both diffusive and bending wave regime) well describe the ing a shadow onto the outer regions. The misaligned inner evolutionofthedisconshorttimescales.Onlonger(closeto disc could be tracing the inclination of the stellar compan- viscous)timescales,itispreferabletomodelthewarpevolu- ion(Billeretal.2012;Closeetal.2014)observedwithinthe tionwith3Dhydrodynamicalsimulations,sincetheyinclude hugemmcavityofthetransitiondisc.Furthermore,CO(6- all the corrected physics in a simple α-disc (Nixon & King 5) ALMA observations of the system have shown a twisted 2016).Amoredetaileddiscussiononsuch3Dsimulationsis velocity centroid map (Casassus et al. 2015), which can be reported in Section 2. modelled by extreme disc warping, probably caused by a Whereas new observations are probing possible signa- dynamicalinteractionwiththelowmasscompanion(Casas- turesofdiscwarping,andthetheoreticalbasisofsuchmech- sus2016).Non-Keplerianflowshavealsobeendetectedwith anism is well established (at least for α-discs), the connec- lower angular resolution in CO (3-2) (Casassus et al. 2013; tionbetweenthetwohasbeenpoorlyaddressed.Inparticu- Rosenfeld et al. 2014). Finally, the KH 15D binary system lar,notmuchattentionhasbeendrawntoobservationaldi- has been modelled as a compact precessing misaligned cir- agnosticsofwarpedprotoplanetarydiscs,inordertoextract cumbinary disc, which causes the long-term modulation of quantitativeinformationwhensuchstructuresareobserved. the light curve of the system (Chiang & Murray-Clay 2004; The effect of the asymmetric illumination of a warped disc Herbst et al. 2010; Lodato & Facchini 2013; Windemuth & was studied by Terquem & Bertout (1996) in circumbinary Herbst 2014). All these observations go together with com- discs while Nixon & Pringle (2010) studied the effect of ir- plementaryrecentfindings,wherediscwarpinghasnotbeen radiation in discs warped due to tidal interactions during detectedyet,butwherecircumprimary(andsometimescir- a stellar fly-by. These studies noticed that warping induces cumsecondary) discs are clearly misaligned with the orbital significantchangesinthespectralenergydistribution(SED) plane of the binary (e.g. Roccatagliata et al. 2011; Jensen atwavelengthsλ>100µm.Similarly,Flaherty&Muzerolle & Akeson 2014; Williams et al. 2014). Some of these discs (2010) provided a prediction for the SEDs of warped discs, areverylikelytoshowkinematicsignaturesofdiscwarping, where the warping was not computed via hydrodynamical when higher angular resolution observations will be avail- simulations, but was simply parametrised in a razor-thin able. Finally, the occurrence of particularly deep dimming model.Rosenfeldetal.(2014)(seealsoCasassusetal.2015) events in the optical of some disc-bearing systems (as AA predicted that warps would cause typical twisted first mo- TauandRWAur)hasbeenalsointerpretedasduetowarp- ment maps in the molecular line emission. Finally, Ruge ing of the inner disc (< few au), where the warp could be etal.(2015)post-processed3Dhydrodynamicalsimulations caused both by magnetic fields (see Bouvier et al. 2013, for of both coplanar and warped circumbinary discs, showing AATau)orunresolvedbinaries(seeFacchinietal.2016,for how the latter cause an asymmetric illumination that can RW Aur). be traced in continuum intensities maps by ALMA. Thetheoryofdiscwarpinghasanticipatedtheobserva- In this paper, we aim to systematically derive observa- tionsbyseveraldecades(intheprotoplanetaryregime).For tionaldiagnosticsofwarpedcircumbinarydiscs.Thewarped arecentreviewonthetheoreticalbasisofthephenomenon, structuresareestimatedusing3Dsmoothedparticlehydro- seeNixon&King(2016).Theresponseofadisctoanonax- dynamics(SPH)simulations,andthenpostprocessedtoob- isymmetric(effective)potentialhasbeeninitiallystudiedby tain both total intensity maps at submm wavelengths, and Papaloizou & Pringle (1983). This study derived the equa- radially resolved molecular emission lines. We also provide tionsfortheevolutionofawarpinadiffusiveregime,which similardiagnosticsintotalintensitymapsinscatteredlight. occurswheneverviscousforcesdominateoverpressureforces We cover a large parameter space, ranging over viscosity, inthedisc.HerewedefinetheaspectratioofadiscasHp/r, inclinations and position angles. These predictions can be whereHpisthescaleheightofthediscandrtheradialcoor- used to interpret present and future observations. dinate,andweparametrisethedisckinematicviscositywith the simple prescription by Shakura & Sunyaev (1973), i.e. ν=αcsHp,wherecs isthesoundspeed,andα adimension- less parameter 1. The criterion for a diffusively evolving (cid:28) disc becomes then α >Hp/r (Papaloizou & Pringle 1983). 2 HYDRODYNAMICAL SIMULATIONS Ogilvie (1999) has extended the analytic approach to a full non-linear regime in the diffusive case. When α<Hp/r, the Inordertoobtainthewarpedstructureofamisalignedcir- warp propagates as a bending wave with a velocity cs/2 cumbinary disc, we use the 3D SPH code phantom1 (Price ∼ (e.g. Papaloizou & Lin 1995; Pringle 1999), since pressure &Federrath2010;Lodato&Price2010;Price2012),which forces dominate the communication of the external torque. has shown significant agreement with both linear and non Linearised equations modelling the warp evolution in such linear theory of warp propagation (Lodato & Price 2010; regimehavebeenderivedbyPapaloizou&Lin(1995),Demi- Facchini et al. 2013; Nealon et al. 2015), and has already anski&Ivanov(1997)andLubow&Ogilvie(2000).Anini- been used to model circumbinary discs (Nixon 2012; Nixon tialstudyonthenon-linearresponseofthesebendingwaves et al. 2013; Facchini et al. 2013; Cazzoletti et al. 2016; Ra- is reported in Ogilvie (2006). For typical physical param- gusa et al. 2017). eters, the diffusive regime describes accretion discs around blackholes,whereasthebendingwavesregimeisapplicable toprotoplanetaryandprotostellardiscs.Manyoftheseana- lyticresultshavebeenconfirmedvianumericalsimulations. In particular, in a linear regime, 3D hydrodynamical sim- 1 http://users.monash.edu.au/~dprice/phantom MNRAS000,1–18(2016) Observability of warps in circumbinary discs 3 2.1 Setup Table 1. Setup of the SPH simulations. All the simulations are For the simulations, we use a setup very similar to the one run for an equal mass binary with separation a, rout=15a, and N=106 particles. reported in Facchini et al. (2013). All the simulations are initiated with two equal mass stars on a circular orbit with a binary separation a. For a coplanar circumbinary disc, α β0 (◦) Binaryorbits N the tidal truncation radius is 1.7a (Artymowicz & Lubow 0.02 15 1130 106 ∼ 1994),anditisnotexpectedtoshrinksignificantlyformod- 0.02 30 1130 106 eratemisalignments(Lubowetal.2015).Wethereforemodel 0.05 15 600 106 the stars with sink particles (e.g. Bate et al. 1995), with an 0.05 30 600 106 accretion radius of 0.95a, and a total binary mass of unity 0.1 15 600 106 0.1 30 600 106 incodeunits.Wearenotinterestedinfollowingthegaseous 0.2 15 600 106 streams flowing onto the binary stars, and we choose such 0.2 30 600 106 large accretion radii to speed up the computation. The ini- tial surface density of the disc is: (cid:18) (cid:114) (cid:19) In the simulations we still have an artificial viscosity, in or- r Σ(r)=Σ0r−1 1 in , (1) der to deal with potential shocks and prevent particle in- − r terpenetrationusingtheMorris&Monaghan(1997)switch. whererin=1.5aandΣ0isanarbitraryscaleparameterinthe In order for the physical viscosity to be higher than the ef- hydrodynamicalsimulations(weneglectselfgravity).Weset fective term generated by the artificial viscosity (Lodato & an initial outer radius rout=15a. The disc is assigned a lo- Price 2010) our simulations have quite high α values. cally isothermal equation of state, with cs∝r−1/4, in order Asmentionedabove,thediscisinitiallyflatonaplane to have a temperature profile that scales as r 1/2. Tem- misaligned to the binary orbit. The disc then evolves, re- − ∼ perature and surface density profiles were chosen to match sponding to the external torque of the central binary, until typical fitted profiles of protoplanetary discs (e.g. Andrews it reaches a warped quasi steady state. The simulations are & Williams 2005, 2007). The aspect ratio Hp/r is set to 0.1 run long enough (for 600 binary orbits for most cases, at r=rin. All the discs are simulated with N =106 SPH 1130 for the less visc∼ous simulations) that both the sur- ∼ particles. face density and the warp reach such quasi steady state. Somesimulations(e.g.Nixonetal.2013;Facchinietal. Thus the warped structure of the disc is maintained un- 2013) have shown that a circumbinary disc can tear apart til the end of our simulation. In Table 1 we summarise the (the same behaviour has been observed both with other mainparametersoftheSPHsimulations.Wedidnotexplore codesandindifferentcontexts,suchascircumprimarydiscs, alargeparameterspacehere,sincewepreferreddoingsoin accretiondiscsaroundspinningblackholes,seee.g.Larwood the post-processing phase for the azimuthal angles and the et al. 1996; Larwood & Papaloizou 1997; Fragner & Nelson inclinations along the line of sight of the simulated discs. 2010; Lodato & Price 2010; Nixon et al. 2012; Dog˘an et al. 2015; Nealon et al. 2015), when the local precession torque exceeds the local viscous torque (in a diffusive regime). In 2.2 Results the bending wave regime, discs seem to tear apart when A warped disc can be parametrised by two angles, which the local precession period induced by the external torque arefunctionofthe(spherical)radiusr(seeFigure1).Inthe is longer than the wave communication timescale (see the thindiscapproximation,wecanassumethatthedisciscom- Appendix in Nixon et al. 2013; Nealon et al. 2015). Even posed by a series of flat, infinitesimally thin rings, each of though the details of tearing in circumbinary discs are yet whichcanbeorientedarbitrarilyinspace.Theorientationof to be thoroughly investigated in the bending wave regime, eachringcanbedescribedbyitsspecificangularmomentum we can confidently prevent the disc from tearing apart, by l(r). To describe the orientation of l(r) in the 3D cartesian simulatingsmall/moderatemisalignmentanglesbetweenthe space we chose a frame of reference in which the angular binary orbit and the plane of the disc. We run simulations momentum vector of the central binary is aligned with the with two different initial misalignment angles, denoted by z-axis. In this case the specific angular momentum of each β , of 15 and 30 . We have checked a posteriori that none 0 ◦ ◦ ringcanbewritteninthefollowingform(e.g.Pringle1996): of our simulations shows any sign of disc breaking. l(r)=(cosγ(r)sinβ(r),sinγ(r)sinβ(r),cosβ(r))T. The tilt an- In order to precisely control the physical viscosity in gle,β(r),definestheanglebetweenthedirectionofthespe- our simulations we use the setup described in section 6.1 of cific angular momentum of the disc and the binary angular Facchinietal.(2013).Wecomputedirectlythestresstensor momentum (i.e. the positive z-axis), while the twist, γ(r), in the Navier-Stokes equation, following the formulation by describestheazimuthalangleofthespecificangularmomen- Flebbeetal.(1994),whichhasprovedtomimicthephysical tum with respect to an arbitrary axis perpendicular to the viscositywithhighprecision(Lodato&Price2010).Weset angularmomentumofthebinary(i.e.inthexy-plane).Thus, the shear viscosity by using the standard α parameter, and themisalignmentbetweenthebinaryorbitandaplanardisc, we set the bulk component to 0. The kinematic viscosity ν as set in our initial conditions of the SPH simulations, can is computed via the equation: be described as a radially independent non-zero tilt angle β .Inordertoobtainβ(r),γ(r)aswellastheradialsurface 0 c2 density profile Σ(r) from the SPH simulations we compute ν=α s (2) azimuthallyaverageddiscquantitiesinthinsphericalshells. Ω MNRAS000,1–18(2016) 4 Juh´asz, Facchini y t Symmetricdisc a -5 nsi e d g o l -6 -7 l z y 10 a -8 x Figure2.Crosssectionofthefinalsnapshotofthemostviscous disc(α=0.2).Colourcodingshowsdensityonalogarithmicscale, Warpeddisc(binaryframe) b inarbitraryunits. Theprocedureisdescribedinsection3.2.6ofLodato&Price (2010). ate InFigure2weshowthecrosssectionofthemostviscous din disc(α=0.2)attheendofthesimulation.Asexpected,the or discshowsawarpintheinnerregions.Tobemorequantita- o c tive,weillustratethetiltandthetwistanglesasafunction l z ofradiusforalloursimulationsinFigure3.Theangleshave z beenestimatedouttor=30a,sincethedischassignificantly β y spreadintheradialdirectionduetopressureforces(theini- γ tial condition truncates the surface density at r=15a) and x viscous spreading. The tilt dependence on radius does not vary much with viscosity, for a given fixed initial misalign- ment(β ).However,thetiltasafunctionofradiusshowsa 0 Warpeddisc c strongdependence( linear)onβ .Thisindicatesthatthe 0 ∼ (frameofdiscouterradius) tiltstructureismostlydeterminedbythenon-Kepleriandisc profile (see Section 3.3 of Foucart & Lai 2013). The semi- analyticalsimulationsbyFacchinietal.(2014)alsoconfirm thatthisshouldhappenintheregionoftheparameterspace explored in this paper. For the highest viscosities, the disc is slowly getting aligned with the binary plane, as expected from viscous damping (Bate et al. 2000; King et al. 2013). On the contrary, the twist angle does not depend on the l z’ theinitialmisalignment,β0.ThisisalsoshowninFigure3b, where we plot the difference between the twist angle in the β0 y’ outer and in the inner regions of the simulated disc (∆γ) as γ0 a function of viscosity. The angle ∆γ is roughly linear with x’ α, as analytically predicted by Foucart & Lai (2013) and - confirmed by Facchini et al. (2014) (in the bending wave regime). Figure 1.Illustrationofthewarpstructurediscussedinthepa- 2.3 Observability of warp fine structures per. The bottom left cartoon shows the orientation of the coor- dinate system while the blue arrow shows the specific angular While the tilt and the twist angles contain important in- momentum vector at the inner edge of the disc. Panel a: The formations on the formation and evolution of warped discs structure of an unperturbed disc whose rotation axis is aligned theyareobservationallyextremelychallengingtodetermine withthez-axisofacartesiancoordinatesystemisaxisymmetric with current instrumentation. The most important reason aroundthez-axisandhasamirrorsymmetrywithrespecttothe forthatisthatβ andγdescribetheorientationofthespecific x-yplane.Panelb:Awarpeddiscischaracterisedbytwoangles: angularmomentumofthediscateachsphericalradiuswith thetwist(β(r))andthetilt(γ(r))definedinacoordinatesystem respect to that of the central binary. Therefore, to observa- alignedwiththeorbitalplaneofthecentralbinary.Thetiltangle tionallyconstrainβ andγ onehastoknowtheorientationof istheanglebetweenthepositivex-axisandtheprojectionofthe thebinaryorbitand/ortheirspecificangularmomentumin specific angular momentum vector onto the x-y plane. Panel c: Without the knowledge of the binary orbit one can only deter- mine quantities defined in a coordinate system aligned with the MNRAS000,1–18(2016) discitselfatacharacteristicradius.Inthisfigurethecoordinate systemisalignedwiththeouteredgeofthedisc. Observability of warps in circumbinary discs 5 30 10 a b 0 25 10 20 − []◦ []◦ 20 β15 ∆γ− 30 − 10 40 − 5 50 − 40 20 c d 30 20 15 10 []◦ []◦ 0 β010 ∆γ0 10 − 5 20 − 30 − 0 40 0 5 10 15 20 25 30 − 0 5 10 15 20 25 30 r/a r/a Figure 3. Tilt and twist angles of the SPH simulations. Colours indicate different viscosities: blue, black, red and green lines are associated with α=0.02, 0.05, 0.1 and 0.2, respectively. Thick lines show the simulations with initial inclination i=30◦, and thin lines simulationswithi=15◦.Panelaandbshowsthetilt(β)andthevariationoftwistbetweentheinnerandtheouterradius(∆γ)defined inthecoordinatesystemalignedtotheorbitofthebinary.Itiseasilyseenthatforafixedbinaryparametersandorbitβ dependsvery weaklyontheviscositybutinsteadisdeterminedbytheinitialmisalignmentanglebetweenthediscandthebinary.Incontrast,∆γ,the totalvariationofγ inthedisc,issetbytheviscosityandisindependentoftheinitialmisalignmentifthebinaryorbitandmassratiois fixed.Panelcanddshowsthetilt(β(cid:48))andthevariationofthetwist(∆γ(cid:48))anglesdefinedinacoordinatesystemalignedtotheplaneof theouteredgeofthedisc(seeSection2.3andFigure1fordefinition).Whilethemostimportantparameterdeterminingβ(cid:48) isstillthe initialmisalignmentthescatterislargerthanforβ inPanela.∆γ(cid:48) seemstodependextremelyweaklyontheviscosity. anabsolutesense.Thisisextremelychallengingtomeasure where represents matrix multiplication, β and γ are out out · as either the binary cannot be spatially resolved, or there the twist and the tilt angles taken at the outer radius of mightbenotenoughastrometricmeasurementstoprecisely the disc, while Ry(βout) and Rz(γout) are rotation matrices determine the orientation of the orbit. around the y- and z-axis, respectively Giventhelackofinformationonthebinaryangularmo-   cosβout 0 sinβout mentum vector, spatially resolved observations of circumbi- − Ry(βout)= 0 1 0  (4) nary discs can only provide us with information in a refer- sinβout 0 cosβout ence frame that is aligned with the disc itself at a certain reference radius. To derive observable quantities in this ref-   cosγout sinγout 0 erenceframeweassumeacoordinatesystemwhosepositive Rz(γout)= sinγout cosγout 0. (5) z-axisisalignedwiththespecificangularmomentumofthe − 0 0 1 outer edge of the disc. We will denote the quantities in this Thus the normalised specific angular momentum com- coordinate system with . Thus β and γ will now describe (cid:48) (cid:48) (cid:48) ponents in the coordinate system aligned with the plane of thewarpstructurewithrespecttothediscitself.β describes (cid:48) the outer edge of the disc will be the relative inclination of a given annulus in the disc with respecttothatoftheouteredgeofthedisc,whileγ defines (cid:48) the direction of the line of nodes of a given annulus with lx(cid:48) = cosβoutcosγoutsinβ(r)cosγ(r)+ respecttoanarbitrarydirectioninthex(cid:48)−y(cid:48) plane.Similar cosβoutsinγoutsinβ(r)sinγ(r) sinβoutcosβ(r) (6) tothedefinitionofβ andγ,β andγ canbecalculatedfrom − the specific angular momentu(cid:48)m at e(cid:48)ach radius. The trans- ly(cid:48) = cosγoutsinβsinγ(r)−sinγoutsinβcosγ(r) (7) formationbetweentheframeofthebinaryandtheframeof lz(cid:48) = sinβoutcosγoutsinβ(r)cosγ(r)+ the outer edge of the disc is given by sinβoutsinγoutsinβ(r)sinγ(r) cosβoutcosβ(r).(8) − Thetwistandthetiltinthisframeofreferencecannow l(cid:48)(x(cid:48),y(cid:48),z(cid:48))=Ry(βout) Rz(γout) l(x,y,z) (3) be calculated from · · MNRAS000,1–18(2016) 6 Juh´asz, Facchini with the plane of an annulus in the disc at a given radius. The transformation between the two frames is given by cosβ(cid:48) = lz(cid:48) (9) tanγ(cid:48) = ly(cid:48). (10) x(cid:48)(cid:48)  cosβ(cid:48) 0 sinβ(cid:48)cosγ(cid:48) sinγ(cid:48) 0x(cid:48) l − x(cid:48) y(cid:48)(cid:48)= 0 1 0 sinγ(cid:48) cosγ(cid:48) 0y(cid:48) We present β and γ as a function of radius in Fig- z(cid:48)(cid:48) sinβ(cid:48) 0 cosβ(cid:48) 0 0 1 z(cid:48) (cid:48) (cid:48) − ure 3c,d. The total tilt of the disc, i.e. the difference in the (12) tilt angle at the outer and the inner radii, is very similar Here,γ andβ aretheradialdependenttwistandtiltangles, bothintheframeofthebinary(∆β)andintheframeofthe (cid:48) (cid:48) respectively, derived from the SPH simulations. The con- disc (∆β(cid:48)). Similar to β the β(cid:48) curves also form two groups nectionbetweenthecartesiancoordinatesandthespherical separatedbytheinitialmisalignmentanglebetweenthecen- coordinates of the spatial grid is given by the usual trans- tralbinary.Thissuggeststhatinoursimulationstheinitial formation misalignment angle is the most important parameter deter- mining the tilt of the disc, and this is independent of the x(cid:48) = rsinθcosφ choice of the frame of reference. Interestingly, ∆γ(cid:48) is signifi- y(cid:48) = rsinθsinφ cantly smaller than ∆γ. The total twist in the frame of the disc is not larger than 10 in any of the models. z(cid:48) = rcosθ. (13) ◦ We assumed that the pressure scale height to be a power-law as a function of radius such that Hp(r) = 0.1(r/r )ζ withr beingtheinnerradiusofthedisc,which in in 3 RADIATIVE TRANSFER weassumedtobeat10au,andζ=0.25istheflaringindex, inagreementwiththesoundspeedradialprofileusedinthe To study the detectability of warps in protoplanetary discs weusedthe3Dradiativetransfercoderadmc-3d2.Thera- SPHsimulations.Thediscouterradiusroutisthusat150au. We assumed that the stars have the parameters of Herbig diative transfer simulations contained two steps. First the stars,i.e.R =2.5R ,M =2M ,T =9500K,andadistance (cid:63) (cid:63) (cid:63) temperatureofthedustinthediscwascalculatedinather- of the source from(cid:12)the observe(cid:12)r of 100pc. Dust particles in mal Monte Carlo simulations, then we calculated images in themodelhadagrainsizedistributionbetween0.1µmand near-infraredscatteredlightaswellasincontinuumandgas 1mm with a power exponent of -3.5. The dust opacity was lines at sub-millimetre wavelengths. calculatedfromtheopticalconstantsofastronomicalsilicate We used a 3D spherical mesh to discretise our model (Weingartner & Draine 2001) using Mie-theory. The disc with N =200, N =180, N =180 grid cells in the radial, r θ φ masswasassumedtobe0.01M withadust-to-gasratioof poloidal and azimuthal directions, respectively. The inner 0.01. The stellar parameters as(cid:12)well as the most important and outer radius of the mesh was chosen to be 10au and parametersofourdiscmodelaresummarizedinTable2.For 150auandthegridcellsweredistributedlogarithmicallyin theCOabundance,withrespecttomolecularhydrogen,we the radial directions. In the poloidal directions 10,80,80,10 assumed 10−4 and a 17O/16O isotopic ratio of 2160. In the grid cells were distributed in an equidistant grid in the [0, wholepaper,withCOwemeanthe12COisotopologue.CO π/2-1.1],[π/2-1.1,π/2],[π/2,π/2+1.1],[π+1.1,π]intervals, anditsisotopologuesareknowntofreezeouttodustgrains respectively. For the azimuthal grid we used an equidistant if the dust temperature drops below 19K decreasing the grid in the [0, 2π] interval. ∼ gas phase abundance of these molecules by several orders Since SPH simulations in general do not have high of magnitude. However, in our models the temperature is enoughresolutionhighabovethemid-plane,wheretheden- above 19K everywhere in the disc due to the high luminos- sity is significantly lower than in the disc mid-plane (thus ity of the stars, thus freeze out of CO and its isotopologues thesmoothinglengthisverylarge),wedidnotusetheden- does not have any effect on our models. To simulate pho- sityoftheSPHsimulationsdirectly.Insteadweusedthetilt todissociation by the UV radiation of the stars we removed (β) and twist (γ) angles as well as the surface density as a all CO from the upper layers of the disc where the radial function of radius (Σ(r)) to generate the density structure opticaldepthasseenfromthestarsat0.5µmislowerthan of the disc in the following way. We described the density unity. structure as We calculate continuum images in the H-band at ρ(x(cid:48)(cid:48),y(cid:48)(cid:48),z(cid:48)(cid:48))= Hp(Σx((cid:48)(cid:48)x,(cid:48)y(cid:48),(cid:48)(cid:48)y)(cid:48)√(cid:48))2πexp(cid:34)−2Hp(zx(cid:48)(cid:48)(cid:48)(cid:48)2,y(cid:48)(cid:48))2(cid:35) (11) 1nst.e6ul5dmyµmatphseaneinffdetc8ht8e0ofJµ=omp3t-i2(cAallLindMeeApotfhBCoaOnndt(h37e4)5li.na7se95pw9r8eol9filGleaHswzec)h.aalTsnoo- calculated channel maps in the J=3-2 transition of 13CO where Hp is the pressure scale height and Σ is the surface (330.5879652218GHz), C18O (329.3305525GHz) and C17O mass density. x, y and z are the cartesian coordinates in a (cid:48) (cid:48) (cid:48) (337.0611298GHz).Forthechannelmapsweassumedave- frameofreferenceinwhichthexy-planeisalignedwiththe (cid:48) (cid:48) locity resolution of 0.2km/s. To study the effect of inclina- planeoftheouteredgeofthedisc.x ,y ,z arethecartesian (cid:48)(cid:48) (cid:48)(cid:48) (cid:48)(cid:48) tionwecalculateimagesat0 ,45 and90 inclinationangles. coordinates in the frame in which the x y -plane is aligned ◦ ◦ ◦ (cid:48)(cid:48) (cid:48)(cid:48) Sincetheinclinationangleofthediscchangesasafunction ofradiusduetothewarpingofthedisc,unlessspecifiedoth- erwise, inclination means the inclination angle at the outer 2 http://www.ita.uni-heidelberg.de/dullemond/software/radmc- edgeofthedisc.Syntheticobservationswerecreatedbycon- 3d/ volvingthenear-infraredimageswitha0.04(cid:48)(cid:48)GaussianPSF, MNRAS000,1–18(2016) Observability of warps in circumbinary discs 7 Table2.Summaryofthestellaranddiscparametersusedinthe 80 14 − radiativetransfercalculations. 60 a 15 − M(cid:63) 2M 60.0 R(cid:63) 2.5R(cid:12) 40 −16 T 9500K(cid:12) MRRMeioddnfufiutssct/(Mdugasst+gas) 0110..05000a11uaMu(cid:12) z[AU]23−83.162202000.704.3 16270.21.9 131.4 95.7 −−−111987logρ Hp(Rin) 0.1 40 20 ζ 0.25 − − 60 21 − − 80 22 − 20 40 60 80 100 − which is the resolution of SPHERE on VLT, the current R[AU] state-of-the-art near-infrared camera. In the sub-millimetre 14 − syntheticobservationswithALMAweregeneratedusingthe b 15 Common Astronomy Software Applications (CASA) v4.2.2 − uansidngthtehcelseiamnotbassekrfvoeritmasakgitnog.gTenheerafutells1y2nmthAetrircayviwsiabsiluitsieeds 0.5 253.3 −16 froersutlhteinsgiminulaatsepdaotibaslerrevsaotliuotnisoninotfwaopdpirffoexriemnattceolynfi0g.u09ra(cid:48)(cid:48)t.iFonosr [rad]θ 0.0 280.0 226.7 280.0 200.0 −1187gρ cuounmtinbuaunmdwsiidmthulaotfiothnes AweLMasAsucmoerrtehlaetofurlaln7d.5cGalHcuzlactoenttihne- /2π− 226.7 200.0 253.3 −−19lo visibilitiesforatotalintegrationtimeof30minutes.Forthe gas line simulations we assumed an integration time of 4h. −0.5 −20 We did not add any instrumental or atmospheric noise 21 − to our synthetic observations. The noise seen in the images 22 arenumericalnoiseoriginatinginhydrodynamicand/orra- 0 1 2 3 4 5 6 − diative transfer calculations or due to the incomplete sam- φ[rad] pling of the Fourier-plane in the synthetic ALMA observa- tions. The reason for the lack of any additional observa- c 17 tional uncertainty in our synthetic observations is that our − intention was to show the amount of asymmetry intrinsic 0.5 18 − to a warped disc without any corruption by an instrument 73.3 110.0 91.7 19 ov(cteioaro.ntgnitd.oahinlietnipotdEaenergasparr,ematenhtteid’cotss)en.raostItnifmtmhwtoehees,epriaehnpssessautrurlretamu.imnmeTgdeehtnneaetorasginslieoovcifehlseneatvhrselaeeelctvmoteeoblarfsyieissnrtuhvicciaadhstn,eiyoowsbnoeosmaebittrseshveeeaorlr--ff /2[rad]πθ−0.0 91.575.0 73.535.0 −−−2210logρ 22 tfahceeebffreicgthstnreeslastaedsymtomtehteryp,rwesheinchcemoafyanwoatrsph,oswucuhpaisnsouurr- −0.5 110.0 −23 − predictions but would be visible in deeper observations. 24 While observables are calculated for all hydrodynamic 0 1 2 3 4 5 6 − modelsinthefollowingweshowtheresultsofthesimulation φ[rad] with α=0.2. To identify the signatures of the warp in the syntheticobservations,wealsocalculatedareferencemodel without a warp. To ensure a meaningful comparison in the Figure 4. Density and temperature structure of a warped cir- referencemodelwetookthesurfacedensityfromourfiducial cumbinarydiscinaverticalsliceintherz-planeatφ=0◦(a)and model with α=0.2, but we set β(r) and γ(r) to zero. in the φ−θ surface at r=12au (b) and at r=100au (c). The colour-scaleimageshowsthelogarithmofthedustdensityinthe disc while the white contours show the dust temperature. The black contour shows the location of the radial optical depth of unityat0.5µm. 4 RESULTS 4.1 Disc structure We present the disc density and temperature structure of from the disc mid-plane, where the density is the highest, our fiducial model in Figure 4a–c. The temperature struc- towards the disc atmosphere, where the density is the low- ture of a warped disc is significantly altered compared to est. In a warped disc this is not necessarily the case. While an unperturbed disc. In a symmetric unperturbed passive close to the inner edge of the disc this inverse correlation irradiateddisctheverticaldensityandtemperatureprofiles still holds (Figure 4b), it breaks further out in the disc. It areinverselycorrelatedsuchthatthetemperatureincreases canbeseeninFigure4cthatat100audistancethevertical MNRAS000,1–18(2016) 8 Juh´asz, Facchini + a b c d e n a pl y x- e v o b a ht g ei H Symmetricdisc Warpeddisc CaseA Warpeddisc CaseB Warpeddisc CaseC − 100 1.5 e f g h 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=0◦ i=0◦ i=0◦ i=0◦ − 10−1 1.5 i j k l 1.0 c] e [arcs 0.5 max offset 00..50 /II,νν c− e D 1.0 − 1.5 i=45◦ i=45◦ i=45◦ i=45◦ − 10−2 1.5 m n o p 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=90◦ i=90◦ i=90◦ i=90◦ − 1.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.5 10−3 − − − − − − − − − − − − RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] Figure 5.Morphologyofwarpeddiscsinnear-infraredscatteredlightobservations.Panelashowsacartoonofasymmetricdiscwhile Panel b–d shows a cartoon depicting the structure and orientation of the warp shown in the synthetic images in each column. The colourscaleshowstheheightabovethemid-planeofthediscwithblueshowingstructuresbelowthemid-planewhileredshowsregions above the mid-plane. The second, third and fourth row from the top show the scattered light images at i=0◦, i=45◦ and i=90◦ inclinationangles,respectively.Thesecond,thirdandfourthcolumnfromtheleftshowsthesyntheticH-bandscatteredlightimagesof ourfiducialwarpeddiscmodel,whiletheleftmostcolumnshowstheimagespredictedimagesofanunperturbed,symmetricdisc.Warped discsshowanazimuthallyasymmetricsurfacebrightnessdistributionatallinclinationangles,evenifseenperfectlyfaceon.Thesurface brightnessasymmetryisduetotheshadowingoftheouterdiscbytheinnerregions.Thepositionangleoftheasymmetryissetbythe orientation/positionangleofthewarp. temperature profile is not coupled to the local density, but density at the inner edge of the disc and the temperature shows correlation (i.e. similar azimuthal modulation) with at larger distances from the star suggests that the primary theverticaldensitystructureintheinnerdiscat12au.This causeofthetemperatureperturbationistheself-shadowing means that at 100au radius, the lowest temperature occurs of the disc by its inner edge where it is curved out of the atdifferentheightsanddensitiesatdifferentazimuthangles. plane of the outer disc. Sinceinapassivediscthedusttemperatureisdeterminedby The requirement for the inner disc to cast a shadow the absorbed stellar radiation, the correlation between the at a given radius r is β(cid:48)(rin)≥Hs(r)/r−Hs(rin)/rin, where MNRAS000,1–18(2016) Observability of warps in circumbinary discs 9 Hs(r) the height of the surface layer above the mid-plane. 4.2.1 Near-infrared images In the unperturbed, reference model Hs(r)/r = 0.23 and Current state-of-the-art NIR cameras (e.g. SPHERE/VLT, Hs(r)/r=0.38 at the inner and outer radius of the disc, re- GPI/Gemini, HiCIAO/SUBARU) can probe regions in the spectively. Thus we would expect to see the effect of non- disc down to several au distances from the central star in axisymmetric self-shadowing if β(cid:48)(rin)≥8.6◦, which is the nearbystar-formingregions.Attheseradiithetemperature case in all but two models (see Figure 3). Even though this inthediscistoolowforsignificantthermalemissionatNIR limitmaychangesomewhatdependingontheopticalprop- wavelengths. Thus NIR observations, dominated by scat- erties of the dust grains in the disc as well as on the flar- tered stellar radiation, probe mainly density variations in ingindex,itisreasonabletoassumethatthecriticalβ (r ) (cid:48) in the upper layers of the disc. above which the non-axisymmetric self-shadowing occurs is on the order of Hp/r. Forface-ondiscs,theNIRimageofanunperturbeddisc Due to the strong three dimensional perturbations one shows a circular symmetric surface brightness distribution expectstwokindsofobservationalsignaturesfromawarped (see Figure 5 e). A warped disc on the other hand shows a circumbinarydisc:surfacebrightnessasymmetriescausedby lopsidedbrightnessdistributionsuchthatthedisccanbedi- the self-shadowing of the disc and kinematical asymmetries videdazimuthallyintoabrightilluminatedandadarkshad- due to the globally non-axisymmetric velocity field. In the owedregioneachhavingapproximatelythesameazimuthal following we will investigate both continuum and gas line extentofabout180◦.Inthefiducialmodelwithα=0.2the observationsandlookforthiskindofobservationalsignature peak-to-peakvariationofthesurfacebrightnessatr=0.5(cid:48)(cid:48) is which can help to identify a warp in an observed disc. afactorof12whileitisafactorof16atr=1(cid:48)(cid:48)forface-onori- Sinceawarpeddiscisagenuinethreedimensionalstruc- entation. It is not only the absolute surface brightness that turethesignaturesdependnotonlyontheinclinationofthe is different in a warped disc compared to a symmetric disc, disc, but also on the azimuth angle. Thus we present syn- butalsotheradialpowerexponentofthesurfacebrightness theticobservablesofwarpeddiscsatthreedifferentorienta- distribution.Inoursymmetriccontrolmodeltheradialsur- tion of the disc depending on the azimuth angle, measured facebrightnessfallsoffasR−2.5,whereRisthedistancefrom in the plane of the outer edge of the disc. The three cases the central star in the plane of the sky at face-on inclina- can be described as follows: tion.Inourfiducialwarpmodeltheradialsurfacebrightness showstrongazimuthalvariation.Fittingapowerlawtothe Case A:Thediscbendsalongthelineofsightsuchthat radial surface brightness profile at the position angle of the theanglebetweenthedirectionoftheobserverandthespe- highest and lowest surface brightness results in a R−3.0 and cificangularmomentumofthediscincreaseswithdecreasing R−2.2 profile, respectively. radius(seeFig.5b).Thismeansthatifthediscisviewedat We also note that the position angle of the low- a non-zero inclination angle the effective inclination of the est/highestsurfacebrightnessdependsontheazimuthangle disc increases with decreasing radius. of the warp. The strong asymmetry in the surface bright- Case B:Thediscbendsalongthelineofsightsuchthat ness as well as its variation with the azimuth angle of the theanglebetweenthedirectionoftheobserverandthespe- warp can be understood in terms of the self-shadowing of cific angular momentum of the disc decreases with decreas- the outer regions of the disc by the inner disc due to the ingradius(seeFig.5c).Ifthediscisviewedatanon-zeroin- warping. In the bright side the curvature of the disc is such clinationangletheeffectiveinclinationofthediscdecreases that it allows stronger irradiation of the disc, i.e. it curves with decreasing radius. Case B can be obtained from the upwards if we are looking at the disc from the top. In the position of Case A by of rotating the disc by 180 around dark side, on the other hand, the situation is the opposite, ◦ the rotation axis of the outer edge of the disc. the outer disc curves downwards into the shadow of the in- CaseC:Thediscbendsperpendiculartothelineofsight ner regions. This explanation is further supported by the (see Fig. 5d). The warping of the disc is such that the pro- fact that the amplitude of the surface brightness asymme- jectedspecificangularmomentumofthediscontotheplane try depends on the strength of the warp, more precisely on of the sky is rotated counter-clockwise with decreasing ra- ∆β(cid:48).Thepeak-to-peakvariationofthesurfacebrightnessin dius. Case C can be obtained from Case A by rotating the the azimuthal direction is only about a factor of 2.5 in case disccounter-clockwiseby90 orfromCaseBbyrotatingthe of the weakest warp (α =0.02 and an initial misalignment ◦ disc clockwise around the rotation axis of the outer edge of of 15 , see Figure A1). For the weakest warp in our hydro- ◦ the disc. dynamic simulations, β(cid:48)(rin)=5.9◦, which is very close to Hs(rin)/rin=0.1. This is a lower value than the criterion we derived in Section 4.1 for β (r ) above which the effect of (cid:48) in non-axisymmetricshelf-shadowingisvisible.Thereasonfor thisisthattheazimuthalasymmetryonthedarksideofthe 4.2 Continuum images discdependsontheself-shadowing,however,ontheopposite Due to the lack of kinematical information, in continuum sitethecurvatureofthediscissuchthattheincidentangle images the observable signatures of disc warps is limited to of the stellar radiation is higher than in the unperturbed surface brightness asymmetries. Variation in the observed disc,whichmakesthatsideofthediscbrighter,makingthe intensity can be caused by the variation of the density, of overall contrast higher than we estimated before from the the temperature or a combination of both. The wavelength self-shadowing alone. of observation and the local dust temperature will deter- For intermediate inclinations (see Figure 5 i-l) the sur- minewhichphysicalparametercausesthesurfacebrightness face brightness shows similar characteristics as for face-on asymmetry. discs, meaning that one side of the disc is brighter than MNRAS000,1–18(2016) 10 Juh´asz, Facchini + a b c d e n a pl y x- e v o b a ht g ei H Symmetricdisc Warpeddisc CaseA Warpeddisc CaseB Warpeddisc CaseC − 1.0 1.5 e f g h 0.9 1.0 c] cse 0.5 0.8 ar [ et 0.0 0.7 s off 0.5 ec− 0.6 D 1.0 − 1.5 i=0◦ i=0◦ i=0◦ i=0◦ 0.5 − 1.5 i j k l 0.4 1.0 c] e [arcs 0.5 0.3 max offset 00..50 /II,νν c− 0.2 e D 1.0 − 1.5 i=45◦ i=45◦ i=45◦ i=45◦ − 1.5 m n o p 0.1 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=90◦ i=90◦ i=90◦ i=90◦ − 1.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.5 − − − − − − − − − − − − RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] Figure 6.SameasFigure5,butthecolorscaleimagesshowsyntheticALMAobservationsat340GHz.Thesynthesisedbeamisshown inwhiteinthebottomleftcornerofeachpanel.Noobviousazimuthalasymmetryisvisibleinthesub-millimetrecontinuumimagesof warped discs. The isophotes of an inclined disc are concentric ellipses, whose position angle and aspect ratio may change as a function ofradius.Thestrengthofthiseffectincreaseswithincreasinginclinationangle. the other. Interestingly in a warped disc the position an- andFigure5l.Thisarciscausedbyphotonsscatteredfrom gle of the dark side is determined by the azimuth angle of the lower side of the disc. In Figure 5j the disc on the near thewarp.Thisisincontrasttotheunperturbeddisc,where side curves downwards, so that the outer regions in the up- the asymmetry is always between the near and far side. In pernearsideofthediscarecompletelyintheshadowofthe an unperturbed disc the far side of the disc looks spatially inner parts. In the lower side of the disc, however, the disc more extended due to the projection, but this can change curves outwards of the shadow of the inner regions making due to anisotropic scattering. For strongly forward peaking it more exposed to direct stellar radiation. scatteringphasefunctionsthenear-sideofthediscmightbe NIR images of edge-on discs show two characteristic brighter. This suggests a caution when determining the in- bright horizontal lanes, caused by scattering of photons in clination and the position angle of a warped disc from NIR the optically thin disc atmosphere (see Figure 5m). These images. Another interesting feature is an arc-like structure two bright arcs are separated by a dark lane in the middle, visible beneath the ellipse of the disc as seen in Figure 5j caused by heavy extinction in the disc mid-plane. Since the MNRAS000,1–18(2016)

ebook img

Observational signatures of linear warps in circumbinary discs PDF

8.6 MB·
by  Attila Juhasz
#journals #arxiv
Save to my drive
Quick download
Download
Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Observational signatures of linear warps in circumbinary discs

MNRAS000,1–18(2016) Preprint11January2017 CompiledusingMNRASLATEXstylefilev3.0 Observational signatures of linear warps in circumbinary discs Attila Juha´sz1(cid:63), Stefano Facchini2† 1Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 2Max-Planck-Institut fu¨r Extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany 7 AcceptedXXX.ReceivedYYY;inoriginalformZZZ 1 0 2 ABSTRACT n In recent years an increasing number of observational studies have hinted at the pres- a ence of warps in protoplanetary discs, however a general comprehensive description J of observational diagnostics of warped discs was missing. We performed a series of 0 3D SPH hydrodynamic simulations and combined them with 3D radiative transfer 1 calculations to study the observability of warps in circumbinary discs, whose plane is misaligned with respect to the orbital plane of the central binary. Our numerical ] P hydrodynamic simulations confirm previous analytical results on the dependence of E the warp structure on the viscosity and the initial misalignment between the binary . andthedisc.Tostudytheobservationalsignaturesofwarpswecalculateimagesinthe h continuumatnear-infraredandsub-millimetrewavelengthsandinthepurerotational p transitionofCOinthesub-millimetre.Warpedcircumbinarydiscsshowsurfacebright- - o nessasymmetryinnear-infraredscatteredlightimagesaswellasinopticallythickgas r lines at sub-millimetre wavelengths. The asymmetry is caused by self-shadowing of t s the disc by the inner warped regions, thus the strength of the asymmetry depends on a the strength of the warp. The projected velocity field, derived from line observations, [ showscharacteristicdeviations,twistsandachangeintheslopeoftherotationcurve, 1 from that of an unperturbed disc. In extreme cases even the direction of rotation v appears to change in the disc inwards of a characteristic radius. The strength of the 1 kinematical signatures of warps decreases with increasing inclination. The strength of 1 all warp signatures decreases with decreasing viscosity. 6 2 Key words: accretion, accretion discs — circumstellar matter — protoplanetary 0 discs — hydrodynamics . 1 0 7 1 1 INTRODUCTION Kotze & Charles 2012, and references therein), and some : v cataclysmic variables (summarised by Olech et al. 2009). Protoplanetarydiscsareusuallythoughtofasaxisymmetric i However, very recently, observational features in a few pro- X structures orbiting a single star. However, it is well known toplanetarydiscshavebeeninterpretedassignaturesofdisc that the geometry of such discs can be more complicated. r warping. In particular, the advent of interferometric obser- a In particular, discs in binary systems can be warped. The vations with the Acatama Large Millimeter Array (ALMA) orbitalplaneofthebinaryandtheplaneoftheouterdiscare hasallowedustostartmeasuringprecisegaskinematicsvia likely to be misaligned due to interaction with other stars molecular line profiles, which is one of the most powerful or to accretion of randomly oriented gas (e.g. Bonnell et al. tools to detect disc warping. For example, Rosenfeld et al. 1992;Bateetal.2010).Suchmisalignmentinducesatorque (2012)haveobservedanexcessintheemissionathighveloc- in the disc, leading to a warped 3D structure. ity in the inner disc of TW Hya, which they interpreted as Untilafewyearsago,welackedobservationalevidence a possible inner warp. Another case is HD 100546, where a of warped protoplanetary discs. In binary systems, tilted warphasbeeninvokedtointerpretthetwistedfirstmoment discs had been inferred only around X-ray binaries show- map of CO (3-2) (Pineda et al. 2014), and even earlier as a ing a long term modulation of the periodic light curve (see possibleexplanationofthespiralstructureobservedinscat- teredlight(Quillen2006)andoftheasymmetriclineprofiles fromsingle-dishobservations(Pani´cetal.2010).Averysig- (cid:63) [email protected] nificant case is HD 142527. Marino et al. (2015) have mod- † [email protected] (cid:13)c 2016TheAuthors 2 Juh´asz, Facchini elledthestronglyasymmetricscatteredlightimageofsucha ulations have shown that the warp evolution equations (in protoplanetarydiscwithahighlymisalignedinnerdisc,cast- both diffusive and bending wave regime) well describe the ing a shadow onto the outer regions. The misaligned inner evolutionofthedisconshorttimescales.Onlonger(closeto disc could be tracing the inclination of the stellar compan- viscous)timescales,itispreferabletomodelthewarpevolu- ion(Billeretal.2012;Closeetal.2014)observedwithinthe tionwith3Dhydrodynamicalsimulations,sincetheyinclude hugemmcavityofthetransitiondisc.Furthermore,CO(6- all the corrected physics in a simple α-disc (Nixon & King 5) ALMA observations of the system have shown a twisted 2016).Amoredetaileddiscussiononsuch3Dsimulationsis velocity centroid map (Casassus et al. 2015), which can be reported in Section 2. modelled by extreme disc warping, probably caused by a Whereas new observations are probing possible signa- dynamicalinteractionwiththelowmasscompanion(Casas- turesofdiscwarping,andthetheoreticalbasisofsuchmech- sus2016).Non-Keplerianflowshavealsobeendetectedwith anism is well established (at least for α-discs), the connec- lower angular resolution in CO (3-2) (Casassus et al. 2013; tionbetweenthetwohasbeenpoorlyaddressed.Inparticu- Rosenfeld et al. 2014). Finally, the KH 15D binary system lar,notmuchattentionhasbeendrawntoobservationaldi- has been modelled as a compact precessing misaligned cir- agnosticsofwarpedprotoplanetarydiscs,inordertoextract cumbinary disc, which causes the long-term modulation of quantitativeinformationwhensuchstructuresareobserved. the light curve of the system (Chiang & Murray-Clay 2004; The effect of the asymmetric illumination of a warped disc Herbst et al. 2010; Lodato & Facchini 2013; Windemuth & was studied by Terquem & Bertout (1996) in circumbinary Herbst 2014). All these observations go together with com- discs while Nixon & Pringle (2010) studied the effect of ir- plementaryrecentfindings,wherediscwarpinghasnotbeen radiation in discs warped due to tidal interactions during detectedyet,butwherecircumprimary(andsometimescir- a stellar fly-by. These studies noticed that warping induces cumsecondary) discs are clearly misaligned with the orbital significantchangesinthespectralenergydistribution(SED) plane of the binary (e.g. Roccatagliata et al. 2011; Jensen atwavelengthsλ>100µm.Similarly,Flaherty&Muzerolle & Akeson 2014; Williams et al. 2014). Some of these discs (2010) provided a prediction for the SEDs of warped discs, areverylikelytoshowkinematicsignaturesofdiscwarping, where the warping was not computed via hydrodynamical when higher angular resolution observations will be avail- simulations, but was simply parametrised in a razor-thin able. Finally, the occurrence of particularly deep dimming model.Rosenfeldetal.(2014)(seealsoCasassusetal.2015) events in the optical of some disc-bearing systems (as AA predicted that warps would cause typical twisted first mo- TauandRWAur)hasbeenalsointerpretedasduetowarp- ment maps in the molecular line emission. Finally, Ruge ing of the inner disc (< few au), where the warp could be etal.(2015)post-processed3Dhydrodynamicalsimulations caused both by magnetic fields (see Bouvier et al. 2013, for of both coplanar and warped circumbinary discs, showing AATau)orunresolvedbinaries(seeFacchinietal.2016,for how the latter cause an asymmetric illumination that can RW Aur). be traced in continuum intensities maps by ALMA. Thetheoryofdiscwarpinghasanticipatedtheobserva- In this paper, we aim to systematically derive observa- tionsbyseveraldecades(intheprotoplanetaryregime).For tionaldiagnosticsofwarpedcircumbinarydiscs.Thewarped arecentreviewonthetheoreticalbasisofthephenomenon, structuresareestimatedusing3Dsmoothedparticlehydro- seeNixon&King(2016).Theresponseofadisctoanonax- dynamics(SPH)simulations,andthenpostprocessedtoob- isymmetric(effective)potentialhasbeeninitiallystudiedby tain both total intensity maps at submm wavelengths, and Papaloizou & Pringle (1983). This study derived the equa- radially resolved molecular emission lines. We also provide tionsfortheevolutionofawarpinadiffusiveregime,which similardiagnosticsintotalintensitymapsinscatteredlight. occurswheneverviscousforcesdominateoverpressureforces We cover a large parameter space, ranging over viscosity, inthedisc.HerewedefinetheaspectratioofadiscasHp/r, inclinations and position angles. These predictions can be whereHpisthescaleheightofthediscandrtheradialcoor- used to interpret present and future observations. dinate,andweparametrisethedisckinematicviscositywith the simple prescription by Shakura & Sunyaev (1973), i.e. ν=αcsHp,wherecs isthesoundspeed,andα adimension- less parameter 1. The criterion for a diffusively evolving (cid:28) disc becomes then α >Hp/r (Papaloizou & Pringle 1983). 2 HYDRODYNAMICAL SIMULATIONS Ogilvie (1999) has extended the analytic approach to a full non-linear regime in the diffusive case. When α<Hp/r, the Inordertoobtainthewarpedstructureofamisalignedcir- warp propagates as a bending wave with a velocity cs/2 cumbinary disc, we use the 3D SPH code phantom1 (Price ∼ (e.g. Papaloizou & Lin 1995; Pringle 1999), since pressure &Federrath2010;Lodato&Price2010;Price2012),which forces dominate the communication of the external torque. has shown significant agreement with both linear and non Linearised equations modelling the warp evolution in such linear theory of warp propagation (Lodato & Price 2010; regimehavebeenderivedbyPapaloizou&Lin(1995),Demi- Facchini et al. 2013; Nealon et al. 2015), and has already anski&Ivanov(1997)andLubow&Ogilvie(2000).Anini- been used to model circumbinary discs (Nixon 2012; Nixon tialstudyonthenon-linearresponseofthesebendingwaves et al. 2013; Facchini et al. 2013; Cazzoletti et al. 2016; Ra- is reported in Ogilvie (2006). For typical physical param- gusa et al. 2017). eters, the diffusive regime describes accretion discs around blackholes,whereasthebendingwavesregimeisapplicable toprotoplanetaryandprotostellardiscs.Manyoftheseana- lyticresultshavebeenconfirmedvianumericalsimulations. In particular, in a linear regime, 3D hydrodynamical sim- 1 http://users.monash.edu.au/~dprice/phantom MNRAS000,1–18(2016) Observability of warps in circumbinary discs 3 2.1 Setup Table 1. Setup of the SPH simulations. All the simulations are For the simulations, we use a setup very similar to the one run for an equal mass binary with separation a, rout=15a, and N=106 particles. reported in Facchini et al. (2013). All the simulations are initiated with two equal mass stars on a circular orbit with a binary separation a. For a coplanar circumbinary disc, α β0 (◦) Binaryorbits N the tidal truncation radius is 1.7a (Artymowicz & Lubow 0.02 15 1130 106 ∼ 1994),anditisnotexpectedtoshrinksignificantlyformod- 0.02 30 1130 106 eratemisalignments(Lubowetal.2015).Wethereforemodel 0.05 15 600 106 the stars with sink particles (e.g. Bate et al. 1995), with an 0.05 30 600 106 accretion radius of 0.95a, and a total binary mass of unity 0.1 15 600 106 0.1 30 600 106 incodeunits.Wearenotinterestedinfollowingthegaseous 0.2 15 600 106 streams flowing onto the binary stars, and we choose such 0.2 30 600 106 large accretion radii to speed up the computation. The ini- tial surface density of the disc is: (cid:18) (cid:114) (cid:19) In the simulations we still have an artificial viscosity, in or- r Σ(r)=Σ0r−1 1 in , (1) der to deal with potential shocks and prevent particle in- − r terpenetrationusingtheMorris&Monaghan(1997)switch. whererin=1.5aandΣ0isanarbitraryscaleparameterinthe In order for the physical viscosity to be higher than the ef- hydrodynamicalsimulations(weneglectselfgravity).Weset fective term generated by the artificial viscosity (Lodato & an initial outer radius rout=15a. The disc is assigned a lo- Price 2010) our simulations have quite high α values. cally isothermal equation of state, with cs∝r−1/4, in order Asmentionedabove,thediscisinitiallyflatonaplane to have a temperature profile that scales as r 1/2. Tem- misaligned to the binary orbit. The disc then evolves, re- − ∼ perature and surface density profiles were chosen to match sponding to the external torque of the central binary, until typical fitted profiles of protoplanetary discs (e.g. Andrews it reaches a warped quasi steady state. The simulations are & Williams 2005, 2007). The aspect ratio Hp/r is set to 0.1 run long enough (for 600 binary orbits for most cases, at r=rin. All the discs are simulated with N =106 SPH 1130 for the less visc∼ous simulations) that both the sur- ∼ particles. face density and the warp reach such quasi steady state. Somesimulations(e.g.Nixonetal.2013;Facchinietal. Thus the warped structure of the disc is maintained un- 2013) have shown that a circumbinary disc can tear apart til the end of our simulation. In Table 1 we summarise the (the same behaviour has been observed both with other mainparametersoftheSPHsimulations.Wedidnotexplore codesandindifferentcontexts,suchascircumprimarydiscs, alargeparameterspacehere,sincewepreferreddoingsoin accretiondiscsaroundspinningblackholes,seee.g.Larwood the post-processing phase for the azimuthal angles and the et al. 1996; Larwood & Papaloizou 1997; Fragner & Nelson inclinations along the line of sight of the simulated discs. 2010; Lodato & Price 2010; Nixon et al. 2012; Dog˘an et al. 2015; Nealon et al. 2015), when the local precession torque exceeds the local viscous torque (in a diffusive regime). In 2.2 Results the bending wave regime, discs seem to tear apart when A warped disc can be parametrised by two angles, which the local precession period induced by the external torque arefunctionofthe(spherical)radiusr(seeFigure1).Inthe is longer than the wave communication timescale (see the thindiscapproximation,wecanassumethatthedisciscom- Appendix in Nixon et al. 2013; Nealon et al. 2015). Even posed by a series of flat, infinitesimally thin rings, each of though the details of tearing in circumbinary discs are yet whichcanbeorientedarbitrarilyinspace.Theorientationof to be thoroughly investigated in the bending wave regime, eachringcanbedescribedbyitsspecificangularmomentum we can confidently prevent the disc from tearing apart, by l(r). To describe the orientation of l(r) in the 3D cartesian simulatingsmall/moderatemisalignmentanglesbetweenthe space we chose a frame of reference in which the angular binary orbit and the plane of the disc. We run simulations momentum vector of the central binary is aligned with the with two different initial misalignment angles, denoted by z-axis. In this case the specific angular momentum of each β , of 15 and 30 . We have checked a posteriori that none 0 ◦ ◦ ringcanbewritteninthefollowingform(e.g.Pringle1996): of our simulations shows any sign of disc breaking. l(r)=(cosγ(r)sinβ(r),sinγ(r)sinβ(r),cosβ(r))T. The tilt an- In order to precisely control the physical viscosity in gle,β(r),definestheanglebetweenthedirectionofthespe- our simulations we use the setup described in section 6.1 of cific angular momentum of the disc and the binary angular Facchinietal.(2013).Wecomputedirectlythestresstensor momentum (i.e. the positive z-axis), while the twist, γ(r), in the Navier-Stokes equation, following the formulation by describestheazimuthalangleofthespecificangularmomen- Flebbeetal.(1994),whichhasprovedtomimicthephysical tum with respect to an arbitrary axis perpendicular to the viscositywithhighprecision(Lodato&Price2010).Weset angularmomentumofthebinary(i.e.inthexy-plane).Thus, the shear viscosity by using the standard α parameter, and themisalignmentbetweenthebinaryorbitandaplanardisc, we set the bulk component to 0. The kinematic viscosity ν as set in our initial conditions of the SPH simulations, can is computed via the equation: be described as a radially independent non-zero tilt angle β .Inordertoobtainβ(r),γ(r)aswellastheradialsurface 0 c2 density profile Σ(r) from the SPH simulations we compute ν=α s (2) azimuthallyaverageddiscquantitiesinthinsphericalshells. Ω MNRAS000,1–18(2016) 4 Juh´asz, Facchini y t Symmetricdisc a -5 nsi e d g o l -6 -7 l z y 10 a -8 x Figure2.Crosssectionofthefinalsnapshotofthemostviscous disc(α=0.2).Colourcodingshowsdensityonalogarithmicscale, Warpeddisc(binaryframe) b inarbitraryunits. Theprocedureisdescribedinsection3.2.6ofLodato&Price (2010). ate InFigure2weshowthecrosssectionofthemostviscous din disc(α=0.2)attheendofthesimulation.Asexpected,the or discshowsawarpintheinnerregions.Tobemorequantita- o c tive,weillustratethetiltandthetwistanglesasafunction l z ofradiusforalloursimulationsinFigure3.Theangleshave z beenestimatedouttor=30a,sincethedischassignificantly β y spreadintheradialdirectionduetopressureforces(theini- γ tial condition truncates the surface density at r=15a) and x viscous spreading. The tilt dependence on radius does not vary much with viscosity, for a given fixed initial misalign- ment(β ).However,thetiltasafunctionofradiusshowsa 0 Warpeddisc c strongdependence( linear)onβ .Thisindicatesthatthe 0 ∼ (frameofdiscouterradius) tiltstructureismostlydeterminedbythenon-Kepleriandisc profile (see Section 3.3 of Foucart & Lai 2013). The semi- analyticalsimulationsbyFacchinietal.(2014)alsoconfirm thatthisshouldhappenintheregionoftheparameterspace explored in this paper. For the highest viscosities, the disc is slowly getting aligned with the binary plane, as expected from viscous damping (Bate et al. 2000; King et al. 2013). On the contrary, the twist angle does not depend on the l z’ theinitialmisalignment,β0.ThisisalsoshowninFigure3b, where we plot the difference between the twist angle in the β0 y’ outer and in the inner regions of the simulated disc (∆γ) as γ0 a function of viscosity. The angle ∆γ is roughly linear with x’ α, as analytically predicted by Foucart & Lai (2013) and - confirmed by Facchini et al. (2014) (in the bending wave regime). Figure 1.Illustrationofthewarpstructurediscussedinthepa- 2.3 Observability of warp fine structures per. The bottom left cartoon shows the orientation of the coor- dinate system while the blue arrow shows the specific angular While the tilt and the twist angles contain important in- momentum vector at the inner edge of the disc. Panel a: The formations on the formation and evolution of warped discs structure of an unperturbed disc whose rotation axis is aligned theyareobservationallyextremelychallengingtodetermine withthez-axisofacartesiancoordinatesystemisaxisymmetric with current instrumentation. The most important reason aroundthez-axisandhasamirrorsymmetrywithrespecttothe forthatisthatβ andγdescribetheorientationofthespecific x-yplane.Panelb:Awarpeddiscischaracterisedbytwoangles: angularmomentumofthediscateachsphericalradiuswith thetwist(β(r))andthetilt(γ(r))definedinacoordinatesystem respect to that of the central binary. Therefore, to observa- alignedwiththeorbitalplaneofthecentralbinary.Thetiltangle tionallyconstrainβ andγ onehastoknowtheorientationof istheanglebetweenthepositivex-axisandtheprojectionofthe thebinaryorbitand/ortheirspecificangularmomentumin specific angular momentum vector onto the x-y plane. Panel c: Without the knowledge of the binary orbit one can only deter- mine quantities defined in a coordinate system aligned with the MNRAS000,1–18(2016) discitselfatacharacteristicradius.Inthisfigurethecoordinate systemisalignedwiththeouteredgeofthedisc. Observability of warps in circumbinary discs 5 30 10 a b 0 25 10 20 − []◦ []◦ 20 β15 ∆γ− 30 − 10 40 − 5 50 − 40 20 c d 30 20 15 10 []◦ []◦ 0 β010 ∆γ0 10 − 5 20 − 30 − 0 40 0 5 10 15 20 25 30 − 0 5 10 15 20 25 30 r/a r/a Figure 3. Tilt and twist angles of the SPH simulations. Colours indicate different viscosities: blue, black, red and green lines are associated with α=0.02, 0.05, 0.1 and 0.2, respectively. Thick lines show the simulations with initial inclination i=30◦, and thin lines simulationswithi=15◦.Panelaandbshowsthetilt(β)andthevariationoftwistbetweentheinnerandtheouterradius(∆γ)defined inthecoordinatesystemalignedtotheorbitofthebinary.Itiseasilyseenthatforafixedbinaryparametersandorbitβ dependsvery weaklyontheviscositybutinsteadisdeterminedbytheinitialmisalignmentanglebetweenthediscandthebinary.Incontrast,∆γ,the totalvariationofγ inthedisc,issetbytheviscosityandisindependentoftheinitialmisalignmentifthebinaryorbitandmassratiois fixed.Panelcanddshowsthetilt(β(cid:48))andthevariationofthetwist(∆γ(cid:48))anglesdefinedinacoordinatesystemalignedtotheplaneof theouteredgeofthedisc(seeSection2.3andFigure1fordefinition).Whilethemostimportantparameterdeterminingβ(cid:48) isstillthe initialmisalignmentthescatterislargerthanforβ inPanela.∆γ(cid:48) seemstodependextremelyweaklyontheviscosity. anabsolutesense.Thisisextremelychallengingtomeasure where represents matrix multiplication, β and γ are out out · as either the binary cannot be spatially resolved, or there the twist and the tilt angles taken at the outer radius of mightbenotenoughastrometricmeasurementstoprecisely the disc, while Ry(βout) and Rz(γout) are rotation matrices determine the orientation of the orbit. around the y- and z-axis, respectively Giventhelackofinformationonthebinaryangularmo-   cosβout 0 sinβout mentum vector, spatially resolved observations of circumbi- − Ry(βout)= 0 1 0  (4) nary discs can only provide us with information in a refer- sinβout 0 cosβout ence frame that is aligned with the disc itself at a certain reference radius. To derive observable quantities in this ref-   cosγout sinγout 0 erenceframeweassumeacoordinatesystemwhosepositive Rz(γout)= sinγout cosγout 0. (5) z-axisisalignedwiththespecificangularmomentumofthe − 0 0 1 outer edge of the disc. We will denote the quantities in this Thus the normalised specific angular momentum com- coordinate system with . Thus β and γ will now describe (cid:48) (cid:48) (cid:48) ponents in the coordinate system aligned with the plane of thewarpstructurewithrespecttothediscitself.β describes (cid:48) the outer edge of the disc will be the relative inclination of a given annulus in the disc with respecttothatoftheouteredgeofthedisc,whileγ defines (cid:48) the direction of the line of nodes of a given annulus with lx(cid:48) = cosβoutcosγoutsinβ(r)cosγ(r)+ respecttoanarbitrarydirectioninthex(cid:48)−y(cid:48) plane.Similar cosβoutsinγoutsinβ(r)sinγ(r) sinβoutcosβ(r) (6) tothedefinitionofβ andγ,β andγ canbecalculatedfrom − the specific angular momentu(cid:48)m at e(cid:48)ach radius. The trans- ly(cid:48) = cosγoutsinβsinγ(r)−sinγoutsinβcosγ(r) (7) formationbetweentheframeofthebinaryandtheframeof lz(cid:48) = sinβoutcosγoutsinβ(r)cosγ(r)+ the outer edge of the disc is given by sinβoutsinγoutsinβ(r)sinγ(r) cosβoutcosβ(r).(8) − Thetwistandthetiltinthisframeofreferencecannow l(cid:48)(x(cid:48),y(cid:48),z(cid:48))=Ry(βout) Rz(γout) l(x,y,z) (3) be calculated from · · MNRAS000,1–18(2016) 6 Juh´asz, Facchini with the plane of an annulus in the disc at a given radius. The transformation between the two frames is given by cosβ(cid:48) = lz(cid:48) (9) tanγ(cid:48) = ly(cid:48). (10) x(cid:48)(cid:48)  cosβ(cid:48) 0 sinβ(cid:48)cosγ(cid:48) sinγ(cid:48) 0x(cid:48) l − x(cid:48) y(cid:48)(cid:48)= 0 1 0 sinγ(cid:48) cosγ(cid:48) 0y(cid:48) We present β and γ as a function of radius in Fig- z(cid:48)(cid:48) sinβ(cid:48) 0 cosβ(cid:48) 0 0 1 z(cid:48) (cid:48) (cid:48) − ure 3c,d. The total tilt of the disc, i.e. the difference in the (12) tilt angle at the outer and the inner radii, is very similar Here,γ andβ aretheradialdependenttwistandtiltangles, bothintheframeofthebinary(∆β)andintheframeofthe (cid:48) (cid:48) respectively, derived from the SPH simulations. The con- disc (∆β(cid:48)). Similar to β the β(cid:48) curves also form two groups nectionbetweenthecartesiancoordinatesandthespherical separatedbytheinitialmisalignmentanglebetweenthecen- coordinates of the spatial grid is given by the usual trans- tralbinary.Thissuggeststhatinoursimulationstheinitial formation misalignment angle is the most important parameter deter- mining the tilt of the disc, and this is independent of the x(cid:48) = rsinθcosφ choice of the frame of reference. Interestingly, ∆γ(cid:48) is signifi- y(cid:48) = rsinθsinφ cantly smaller than ∆γ. The total twist in the frame of the disc is not larger than 10 in any of the models. z(cid:48) = rcosθ. (13) ◦ We assumed that the pressure scale height to be a power-law as a function of radius such that Hp(r) = 0.1(r/r )ζ withr beingtheinnerradiusofthedisc,which in in 3 RADIATIVE TRANSFER weassumedtobeat10au,andζ=0.25istheflaringindex, inagreementwiththesoundspeedradialprofileusedinthe To study the detectability of warps in protoplanetary discs weusedthe3Dradiativetransfercoderadmc-3d2.Thera- SPHsimulations.Thediscouterradiusroutisthusat150au. We assumed that the stars have the parameters of Herbig diative transfer simulations contained two steps. First the stars,i.e.R =2.5R ,M =2M ,T =9500K,andadistance (cid:63) (cid:63) (cid:63) temperatureofthedustinthediscwascalculatedinather- of the source from(cid:12)the observe(cid:12)r of 100pc. Dust particles in mal Monte Carlo simulations, then we calculated images in themodelhadagrainsizedistributionbetween0.1µmand near-infraredscatteredlightaswellasincontinuumandgas 1mm with a power exponent of -3.5. The dust opacity was lines at sub-millimetre wavelengths. calculatedfromtheopticalconstantsofastronomicalsilicate We used a 3D spherical mesh to discretise our model (Weingartner & Draine 2001) using Mie-theory. The disc with N =200, N =180, N =180 grid cells in the radial, r θ φ masswasassumedtobe0.01M withadust-to-gasratioof poloidal and azimuthal directions, respectively. The inner 0.01. The stellar parameters as(cid:12)well as the most important and outer radius of the mesh was chosen to be 10au and parametersofourdiscmodelaresummarizedinTable2.For 150auandthegridcellsweredistributedlogarithmicallyin theCOabundance,withrespecttomolecularhydrogen,we the radial directions. In the poloidal directions 10,80,80,10 assumed 10−4 and a 17O/16O isotopic ratio of 2160. In the grid cells were distributed in an equidistant grid in the [0, wholepaper,withCOwemeanthe12COisotopologue.CO π/2-1.1],[π/2-1.1,π/2],[π/2,π/2+1.1],[π+1.1,π]intervals, anditsisotopologuesareknowntofreezeouttodustgrains respectively. For the azimuthal grid we used an equidistant if the dust temperature drops below 19K decreasing the grid in the [0, 2π] interval. ∼ gas phase abundance of these molecules by several orders Since SPH simulations in general do not have high of magnitude. However, in our models the temperature is enoughresolutionhighabovethemid-plane,wheretheden- above 19K everywhere in the disc due to the high luminos- sity is significantly lower than in the disc mid-plane (thus ity of the stars, thus freeze out of CO and its isotopologues thesmoothinglengthisverylarge),wedidnotusetheden- does not have any effect on our models. To simulate pho- sityoftheSPHsimulationsdirectly.Insteadweusedthetilt todissociation by the UV radiation of the stars we removed (β) and twist (γ) angles as well as the surface density as a all CO from the upper layers of the disc where the radial function of radius (Σ(r)) to generate the density structure opticaldepthasseenfromthestarsat0.5µmislowerthan of the disc in the following way. We described the density unity. structure as We calculate continuum images in the H-band at ρ(x(cid:48)(cid:48),y(cid:48)(cid:48),z(cid:48)(cid:48))= Hp(Σx((cid:48)(cid:48)x,(cid:48)y(cid:48),(cid:48)(cid:48)y)(cid:48)√(cid:48))2πexp(cid:34)−2Hp(zx(cid:48)(cid:48)(cid:48)(cid:48)2,y(cid:48)(cid:48))2(cid:35) (11) 1nst.e6ul5dmyµmatphseaneinffdetc8ht8e0ofJµ=omp3t-i2(cAallLindMeeApotfhBCoaOnndt(h37e4)5li.na7se95pw9r8eol9filGleaHswzec)h.aalTsnoo- calculated channel maps in the J=3-2 transition of 13CO where Hp is the pressure scale height and Σ is the surface (330.5879652218GHz), C18O (329.3305525GHz) and C17O mass density. x, y and z are the cartesian coordinates in a (cid:48) (cid:48) (cid:48) (337.0611298GHz).Forthechannelmapsweassumedave- frameofreferenceinwhichthexy-planeisalignedwiththe (cid:48) (cid:48) locity resolution of 0.2km/s. To study the effect of inclina- planeoftheouteredgeofthedisc.x ,y ,z arethecartesian (cid:48)(cid:48) (cid:48)(cid:48) (cid:48)(cid:48) tionwecalculateimagesat0 ,45 and90 inclinationangles. coordinates in the frame in which the x y -plane is aligned ◦ ◦ ◦ (cid:48)(cid:48) (cid:48)(cid:48) Sincetheinclinationangleofthediscchangesasafunction ofradiusduetothewarpingofthedisc,unlessspecifiedoth- erwise, inclination means the inclination angle at the outer 2 http://www.ita.uni-heidelberg.de/dullemond/software/radmc- edgeofthedisc.Syntheticobservationswerecreatedbycon- 3d/ volvingthenear-infraredimageswitha0.04(cid:48)(cid:48)GaussianPSF, MNRAS000,1–18(2016) Observability of warps in circumbinary discs 7 Table2.Summaryofthestellaranddiscparametersusedinthe 80 14 − radiativetransfercalculations. 60 a 15 − M(cid:63) 2M 60.0 R(cid:63) 2.5R(cid:12) 40 −16 T 9500K(cid:12) MRRMeioddnfufiutssct/(Mdugasst+gas) 0110..05000a11uaMu(cid:12) z[AU]23−83.162202000.704.3 16270.21.9 131.4 95.7 −−−111987logρ Hp(Rin) 0.1 40 20 ζ 0.25 − − 60 21 − − 80 22 − 20 40 60 80 100 − which is the resolution of SPHERE on VLT, the current R[AU] state-of-the-art near-infrared camera. In the sub-millimetre 14 − syntheticobservationswithALMAweregeneratedusingthe b 15 Common Astronomy Software Applications (CASA) v4.2.2 − uansidngthtehcelseiamnotbassekrfvoeritmasakgitnog.gTenheerafutells1y2nmthAetrircayviwsiabsiluitsieeds 0.5 253.3 −16 froersutlhteinsgiminulaatsepdaotibaslerrevsaotliuotnisoninotfwaopdpirffoexriemnattceolynfi0g.u09ra(cid:48)(cid:48)t.iFonosr [rad]θ 0.0 280.0 226.7 280.0 200.0 −1187gρ cuounmtinbuaunmdwsiidmthulaotfiothnes AweLMasAsucmoerrtehlaetofurlaln7d.5cGalHcuzlactoenttihne- /2π− 226.7 200.0 253.3 −−19lo visibilitiesforatotalintegrationtimeof30minutes.Forthe gas line simulations we assumed an integration time of 4h. −0.5 −20 We did not add any instrumental or atmospheric noise 21 − to our synthetic observations. The noise seen in the images 22 arenumericalnoiseoriginatinginhydrodynamicand/orra- 0 1 2 3 4 5 6 − diative transfer calculations or due to the incomplete sam- φ[rad] pling of the Fourier-plane in the synthetic ALMA observa- tions. The reason for the lack of any additional observa- c 17 tional uncertainty in our synthetic observations is that our − intention was to show the amount of asymmetry intrinsic 0.5 18 − to a warped disc without any corruption by an instrument 73.3 110.0 91.7 19 ov(cteioaro.ntgnitd.oahinlietnipotdEaenergasparr,ematenhtteid’cotss)en.raostItnifmtmhwtoehees,epriaehnpssessautrurlretamu.imnmeTgdeehtnneaetorasginslieoovcifehlseneatvhrselaeeelctvmoteeoblarfsyieissnrtuhvicciaadhstn,eiyoowsbnoeosmaebittrseshveeeaorlr--ff /2[rad]πθ−0.0 91.575.0 73.535.0 −−−2210logρ 22 tfahceeebffreicgthstnreeslastaedsymtomtehteryp,rwesheinchcemoafyanwoatrsph,oswucuhpaisnsouurr- −0.5 110.0 −23 − predictions but would be visible in deeper observations. 24 While observables are calculated for all hydrodynamic 0 1 2 3 4 5 6 − modelsinthefollowingweshowtheresultsofthesimulation φ[rad] with α=0.2. To identify the signatures of the warp in the syntheticobservations,wealsocalculatedareferencemodel without a warp. To ensure a meaningful comparison in the Figure 4. Density and temperature structure of a warped cir- referencemodelwetookthesurfacedensityfromourfiducial cumbinarydiscinaverticalsliceintherz-planeatφ=0◦(a)and model with α=0.2, but we set β(r) and γ(r) to zero. in the φ−θ surface at r=12au (b) and at r=100au (c). The colour-scaleimageshowsthelogarithmofthedustdensityinthe disc while the white contours show the dust temperature. The black contour shows the location of the radial optical depth of unityat0.5µm. 4 RESULTS 4.1 Disc structure We present the disc density and temperature structure of from the disc mid-plane, where the density is the highest, our fiducial model in Figure 4a–c. The temperature struc- towards the disc atmosphere, where the density is the low- ture of a warped disc is significantly altered compared to est. In a warped disc this is not necessarily the case. While an unperturbed disc. In a symmetric unperturbed passive close to the inner edge of the disc this inverse correlation irradiateddisctheverticaldensityandtemperatureprofiles still holds (Figure 4b), it breaks further out in the disc. It areinverselycorrelatedsuchthatthetemperatureincreases canbeseeninFigure4cthatat100audistancethevertical MNRAS000,1–18(2016) 8 Juh´asz, Facchini + a b c d e n a pl y x- e v o b a ht g ei H Symmetricdisc Warpeddisc CaseA Warpeddisc CaseB Warpeddisc CaseC − 100 1.5 e f g h 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=0◦ i=0◦ i=0◦ i=0◦ − 10−1 1.5 i j k l 1.0 c] e [arcs 0.5 max offset 00..50 /II,νν c− e D 1.0 − 1.5 i=45◦ i=45◦ i=45◦ i=45◦ − 10−2 1.5 m n o p 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=90◦ i=90◦ i=90◦ i=90◦ − 1.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.5 10−3 − − − − − − − − − − − − RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] Figure 5.Morphologyofwarpeddiscsinnear-infraredscatteredlightobservations.Panelashowsacartoonofasymmetricdiscwhile Panel b–d shows a cartoon depicting the structure and orientation of the warp shown in the synthetic images in each column. The colourscaleshowstheheightabovethemid-planeofthediscwithblueshowingstructuresbelowthemid-planewhileredshowsregions above the mid-plane. The second, third and fourth row from the top show the scattered light images at i=0◦, i=45◦ and i=90◦ inclinationangles,respectively.Thesecond,thirdandfourthcolumnfromtheleftshowsthesyntheticH-bandscatteredlightimagesof ourfiducialwarpeddiscmodel,whiletheleftmostcolumnshowstheimagespredictedimagesofanunperturbed,symmetricdisc.Warped discsshowanazimuthallyasymmetricsurfacebrightnessdistributionatallinclinationangles,evenifseenperfectlyfaceon.Thesurface brightnessasymmetryisduetotheshadowingoftheouterdiscbytheinnerregions.Thepositionangleoftheasymmetryissetbythe orientation/positionangleofthewarp. temperature profile is not coupled to the local density, but density at the inner edge of the disc and the temperature shows correlation (i.e. similar azimuthal modulation) with at larger distances from the star suggests that the primary theverticaldensitystructureintheinnerdiscat12au.This causeofthetemperatureperturbationistheself-shadowing means that at 100au radius, the lowest temperature occurs of the disc by its inner edge where it is curved out of the atdifferentheightsanddensitiesatdifferentazimuthangles. plane of the outer disc. Sinceinapassivediscthedusttemperatureisdeterminedby The requirement for the inner disc to cast a shadow the absorbed stellar radiation, the correlation between the at a given radius r is β(cid:48)(rin)≥Hs(r)/r−Hs(rin)/rin, where MNRAS000,1–18(2016) Observability of warps in circumbinary discs 9 Hs(r) the height of the surface layer above the mid-plane. 4.2.1 Near-infrared images In the unperturbed, reference model Hs(r)/r = 0.23 and Current state-of-the-art NIR cameras (e.g. SPHERE/VLT, Hs(r)/r=0.38 at the inner and outer radius of the disc, re- GPI/Gemini, HiCIAO/SUBARU) can probe regions in the spectively. Thus we would expect to see the effect of non- disc down to several au distances from the central star in axisymmetric self-shadowing if β(cid:48)(rin)≥8.6◦, which is the nearbystar-formingregions.Attheseradiithetemperature case in all but two models (see Figure 3). Even though this inthediscistoolowforsignificantthermalemissionatNIR limitmaychangesomewhatdependingontheopticalprop- wavelengths. Thus NIR observations, dominated by scat- erties of the dust grains in the disc as well as on the flar- tered stellar radiation, probe mainly density variations in ingindex,itisreasonabletoassumethatthecriticalβ (r ) (cid:48) in the upper layers of the disc. above which the non-axisymmetric self-shadowing occurs is on the order of Hp/r. Forface-ondiscs,theNIRimageofanunperturbeddisc Due to the strong three dimensional perturbations one shows a circular symmetric surface brightness distribution expectstwokindsofobservationalsignaturesfromawarped (see Figure 5 e). A warped disc on the other hand shows a circumbinarydisc:surfacebrightnessasymmetriescausedby lopsidedbrightnessdistributionsuchthatthedisccanbedi- the self-shadowing of the disc and kinematical asymmetries videdazimuthallyintoabrightilluminatedandadarkshad- due to the globally non-axisymmetric velocity field. In the owedregioneachhavingapproximatelythesameazimuthal following we will investigate both continuum and gas line extentofabout180◦.Inthefiducialmodelwithα=0.2the observationsandlookforthiskindofobservationalsignature peak-to-peakvariationofthesurfacebrightnessatr=0.5(cid:48)(cid:48) is which can help to identify a warp in an observed disc. afactorof12whileitisafactorof16atr=1(cid:48)(cid:48)forface-onori- Sinceawarpeddiscisagenuinethreedimensionalstruc- entation. It is not only the absolute surface brightness that turethesignaturesdependnotonlyontheinclinationofthe is different in a warped disc compared to a symmetric disc, disc, but also on the azimuth angle. Thus we present syn- butalsotheradialpowerexponentofthesurfacebrightness theticobservablesofwarpeddiscsatthreedifferentorienta- distribution.Inoursymmetriccontrolmodeltheradialsur- tion of the disc depending on the azimuth angle, measured facebrightnessfallsoffasR−2.5,whereRisthedistancefrom in the plane of the outer edge of the disc. The three cases the central star in the plane of the sky at face-on inclina- can be described as follows: tion.Inourfiducialwarpmodeltheradialsurfacebrightness showstrongazimuthalvariation.Fittingapowerlawtothe Case A:Thediscbendsalongthelineofsightsuchthat radial surface brightness profile at the position angle of the theanglebetweenthedirectionoftheobserverandthespe- highest and lowest surface brightness results in a R−3.0 and cificangularmomentumofthediscincreaseswithdecreasing R−2.2 profile, respectively. radius(seeFig.5b).Thismeansthatifthediscisviewedat We also note that the position angle of the low- a non-zero inclination angle the effective inclination of the est/highestsurfacebrightnessdependsontheazimuthangle disc increases with decreasing radius. of the warp. The strong asymmetry in the surface bright- Case B:Thediscbendsalongthelineofsightsuchthat ness as well as its variation with the azimuth angle of the theanglebetweenthedirectionoftheobserverandthespe- warp can be understood in terms of the self-shadowing of cific angular momentum of the disc decreases with decreas- the outer regions of the disc by the inner disc due to the ingradius(seeFig.5c).Ifthediscisviewedatanon-zeroin- warping. In the bright side the curvature of the disc is such clinationangletheeffectiveinclinationofthediscdecreases that it allows stronger irradiation of the disc, i.e. it curves with decreasing radius. Case B can be obtained from the upwards if we are looking at the disc from the top. In the position of Case A by of rotating the disc by 180 around dark side, on the other hand, the situation is the opposite, ◦ the rotation axis of the outer edge of the disc. the outer disc curves downwards into the shadow of the in- CaseC:Thediscbendsperpendiculartothelineofsight ner regions. This explanation is further supported by the (see Fig. 5d). The warping of the disc is such that the pro- fact that the amplitude of the surface brightness asymme- jectedspecificangularmomentumofthediscontotheplane try depends on the strength of the warp, more precisely on of the sky is rotated counter-clockwise with decreasing ra- ∆β(cid:48).Thepeak-to-peakvariationofthesurfacebrightnessin dius. Case C can be obtained from Case A by rotating the the azimuthal direction is only about a factor of 2.5 in case disccounter-clockwiseby90 orfromCaseBbyrotatingthe of the weakest warp (α =0.02 and an initial misalignment ◦ disc clockwise around the rotation axis of the outer edge of of 15 , see Figure A1). For the weakest warp in our hydro- ◦ the disc. dynamic simulations, β(cid:48)(rin)=5.9◦, which is very close to Hs(rin)/rin=0.1. This is a lower value than the criterion we derived in Section 4.1 for β (r ) above which the effect of (cid:48) in non-axisymmetricshelf-shadowingisvisible.Thereasonfor thisisthattheazimuthalasymmetryonthedarksideofthe 4.2 Continuum images discdependsontheself-shadowing,however,ontheopposite Due to the lack of kinematical information, in continuum sitethecurvatureofthediscissuchthattheincidentangle images the observable signatures of disc warps is limited to of the stellar radiation is higher than in the unperturbed surface brightness asymmetries. Variation in the observed disc,whichmakesthatsideofthediscbrighter,makingthe intensity can be caused by the variation of the density, of overall contrast higher than we estimated before from the the temperature or a combination of both. The wavelength self-shadowing alone. of observation and the local dust temperature will deter- For intermediate inclinations (see Figure 5 i-l) the sur- minewhichphysicalparametercausesthesurfacebrightness face brightness shows similar characteristics as for face-on asymmetry. discs, meaning that one side of the disc is brighter than MNRAS000,1–18(2016) 10 Juh´asz, Facchini + a b c d e n a pl y x- e v o b a ht g ei H Symmetricdisc Warpeddisc CaseA Warpeddisc CaseB Warpeddisc CaseC − 1.0 1.5 e f g h 0.9 1.0 c] cse 0.5 0.8 ar [ et 0.0 0.7 s off 0.5 ec− 0.6 D 1.0 − 1.5 i=0◦ i=0◦ i=0◦ i=0◦ 0.5 − 1.5 i j k l 0.4 1.0 c] e [arcs 0.5 0.3 max offset 00..50 /II,νν c− 0.2 e D 1.0 − 1.5 i=45◦ i=45◦ i=45◦ i=45◦ − 1.5 m n o p 0.1 1.0 c] e cs 0.5 ar [ et 0.0 s off 0.5 c− e D 1.0 − 1.5 i=90◦ i=90◦ i=90◦ i=90◦ − 1.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.51.5 1.0 0.5 0.0 0.5 1.0 1.5 − − − − − − − − − − − − RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] RAoffset[arcsec] Figure 6.SameasFigure5,butthecolorscaleimagesshowsyntheticALMAobservationsat340GHz.Thesynthesisedbeamisshown inwhiteinthebottomleftcornerofeachpanel.Noobviousazimuthalasymmetryisvisibleinthesub-millimetrecontinuumimagesof warped discs. The isophotes of an inclined disc are concentric ellipses, whose position angle and aspect ratio may change as a function ofradius.Thestrengthofthiseffectincreaseswithincreasinginclinationangle. the other. Interestingly in a warped disc the position an- andFigure5l.Thisarciscausedbyphotonsscatteredfrom gle of the dark side is determined by the azimuth angle of the lower side of the disc. In Figure 5j the disc on the near thewarp.Thisisincontrasttotheunperturbeddisc,where side curves downwards, so that the outer regions in the up- the asymmetry is always between the near and far side. In pernearsideofthediscarecompletelyintheshadowofthe an unperturbed disc the far side of the disc looks spatially inner parts. In the lower side of the disc, however, the disc more extended due to the projection, but this can change curves outwards of the shadow of the inner regions making due to anisotropic scattering. For strongly forward peaking it more exposed to direct stellar radiation. scatteringphasefunctionsthenear-sideofthediscmightbe NIR images of edge-on discs show two characteristic brighter. This suggests a caution when determining the in- bright horizontal lanes, caused by scattering of photons in clination and the position angle of a warped disc from NIR the optically thin disc atmosphere (see Figure 5m). These images. Another interesting feature is an arc-like structure two bright arcs are separated by a dark lane in the middle, visible beneath the ellipse of the disc as seen in Figure 5j caused by heavy extinction in the disc mid-plane. Since the MNRAS000,1–18(2016)

See more

The list of books you might like

Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.
DMCA & Copyright: Dear all, most of the website is community built, users are uploading hundred of books everyday, which makes really hard for us to identify copyrighted material, please contact us if you want any material removed.
Copyright @ PDFdrive.to, 2022 | We love our users