Observational signatures of linear warps in circumbinary discs

8.6 MB·

by  Attila Juhasz

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)