Cyclic Variability of the Circumstellar Disc of the Be Star $ζ$ Tau. II. Testing the 2D Global Disc Oscillation Model PDF
Preview Cyclic Variability of the Circumstellar Disc of the Be Star $ζ$ Tau. II. Testing the 2D Global Disc Oscillation Model
Astronomy&Astrophysicsmanuscriptno.aa6 c ESO2009 (cid:13) January26,2009 Cyclic Variability of the Circumstellar Disc of the Be Star ζ Tau II. Testing the 2D Global Disc Oscillation Model A.C.Carciofi1,A.T.Okazaki2,J.-B.leBouquin3,S.Sˇtefl3,Th.Rivinius3,D.Baade4,J.E.Bjorkman5,andC.A. Hummel4 9 0 0 1 InstitutodeAstronomia,Geof´ısicaeCieˆnciasAtmosfe´ricas,UniversidadedeSa˜oPaulo,RuadoMata˜o1226,CidadeUniversita´ria, 2 Sa˜oPaulo,SP05508-900,Brazil e-mail:[email protected] n 2 FacultyofEngineering,Hokkai-GakuenUniversity,Toyohira-ku,Sapporo062-8605,Japan a J 3 ESO-EuropeanOrganisationforAstronomicalResearchintheSouthernHemisphere,Casilla19001,Santiago19,Chile 8 4 ESO-EuropeanOrganisationforAstronomicalResearchintheSouthernHemisphere,Karl-Schwarzschild-Str.2,85748Garching ] beiMu¨nchen,Germany R 5 UniversityofToledo,DepartmentofPhysics&Astronomy,MS1112801W.BancroftStreetToledo,OH43606USA S . h Received:<date>;accepted:<date>;LATEXed:January26,2009 p - o ABSTRACT r t Context.About 2/3 oftheBestarspresent thesocalledV/Rvariations, aphenomenon characterizedbythequasi-cyclicvariation s a oftheratiobetweenthevioletandredemissionpeaksoftheHIemissionlines.Thesevariationsaregenerallyexplainedbyglobal [ oscillationsinthecircumstellardiscformingaone-armedspiraldensitypatternthatprecessesaroundthestarwithaperiodofafew years. 1 Aims.Inthispaperwemodel,inaself-consistentway,polarimetric,photometric,spectrophotometricandinterferometricobservations v oftheclassicalBestarζTauri.Ourprimarygoalistoconductacriticalquantitativetestoftheglobaloscillationscenario. 8 Methods. We have carried out detailed three-dimensional, NLTE radiative transfer calculations using the radiative transfer code 9 HDUST.Fortheinputforthecodewehaveusedthemostup-to-dateresearchonBestarstoincludeaphysicallyrealisticdescription 0 forthecentralstarandthecircumstellardisc.Weadoptarotationallydeformed,gravitydarkenedcentralstar,surroundedbyadisc 1 whoseunperturbedstateisgivenbyasteady-stateviscousdecretiondiscmodel.Wefurtherassumethatdiscisinverticalhydrostatic . equilibrium. 1 Results. Byadoptingaviscousdecretiondiscmodelforζ Tauriandarigoroussolutionoftheradiativetransfer,wehaveobtained 0 averygoodfitofthetime-averagepropertiesofthedisc.Thisprovidesstrongtheoreticalevidence thattheviscousdecretiondisc 9 modelisthemechanismresponsiblefordiscformation.Withtheglobaloscillationmodelwehavesuccessfullyfittedspatiallyresolved 0 VLTI/AMBERobservationsandthetemporalV/RvariationsoftheHαandBrγlines.Thisresultconvincinglydemonstratesthatthe : v oscillationpatterninthediscisaone-armedspiral.Possiblemodelshortcomings,aswellassuggestionsforfutureimprovements,are i alsodiscussed. X Keywords.Methods:numerical–RadiativeTransfer–Stars:emission-line,Be–Stars:individual:ζTau–Polarization–Techniques: r a interferometric 1. Introduction I, a compilation of spectroscopic, spectrophotometric, polari- metric and spectropolarimetric observations of ζ Tau was pre- ζ Tauriisanorthern,nearbyBestarthathasdrawntheattention sented,covering3completeV/Rcyclesfrom1992tothepresent. ofgenerationsofobserversandtheoreticians.Itspositiononthe Thedatawascomplementedbythefirstspectro-interferometric skymakesitaconvenienttargetforbothnorthernandsouthern observationsofthisstar,madeinDecember,2006withAMBER telescopes,and,asaresult,thisstarhasaveryrichobservational at VLTI. The reader is referred to Paper I for details about the history.ζ Tau hassome distinctivefeaturesthatmake itan im- observationsandinitialdataanalysis. portant laboratory for testing theoretical ideas about processes In this paper we report a detailed modeling of the avail- associatedwiththeBephenomenon.First,ithasacircumstellar able data. We adopt a physical model for the central star and (CS)discthathasshownlittleornosecularevolutioninthepast theCSdisc thatisconsistentwiththemostup-to-dateresearch 18yearsorso(Sˇtefletal.2009,hereafterPaperI;Riviniusetal. onBe stars,asdescribedbelow.Forthecentralstar weadopta 2006).Second,itshowsaverystableV/R1 cycle,withaquasi- model that includes rotational effects such as geometrical de- periodofabout1430days(PaperI; Sˇtefletal.,2007). InPaper formation (DomicianodeSouzaetal., 2003) and the latitudi- Sendoffprintrequeststo:A.C.Carciofi,e-mail:[email protected] 1 TheV/Rtermreferstothequasi-periodicvariationoftheratiobe- It isobserved inζ Tau and numerous other Bestars. SeePaper Iand tweenthevioletandredemissionpeaksofcircumstellaremissionlines. referencesthereinformoredetails. 2 Carciofietal.:Testingthe2DGlobalDiscOscillationModel naldependencyofthe stellar radiation(Townsendetal., 2004). 2. ModelDescription Rotationaleffectsareveryimportantindeterminingthelocalra- ThecomputercodeHDUSTisafully3D,non-localthermody- diation field at a given point in the stellar environment, since namic equilibrium (NLTE) code designed to solve the coupled they redirect part of the flux toward the polar regions. For the problems of radiative transfer, radiative equilibrium and statis- unperturbedstateoftheCSdiscweadoptasteady-state,viscous ticalequilibriumforarbitrarygasdensityandvelocitydistribu- decretion disc model (Leeetal., 1991; Porter, 1999; Okazaki, tions. TheNLTE MonteCarlo simulationperformsa fullspec- 2001;Carciofietal.,2006).Inthismodel,thematerialisejected tralsynthesisbyemittingphotonsfromarotationallydeformed fromthestellarequatorialsurface,driftsoutwardsowingtovis- andgravitydarkenedstar.Thestarisdividedinanumberoflat- cous effects, and forms a geometrically thin disc with nearly itude bins (typically 100), each with its effective temperature, Keplerianrotation(Riviniusetal.,2006). gravity and a spectral shape given by the appropriate Kurucz In order to model both V/R variations of the spectral lines model atmosphere (Kurucz, 1994). After emission by the star, andtheVLTI/AMBERobservations,weusetheglobaldiscos- eachphotonisfollowedasittravelsthroughtheenvelope(where cillationmodelofOkazaki(1997)andPapaloizouetal.(1992). itmaybescattered,orabsorbedandreemitted,manytimes)un- Inthismodel,aone-armedm=1oscillationmodeissuperposed til it escapes. Duringa simulation,whenevera photonscatters, ontheunperturbedstate ofthedisc.Themodelisbasedonthe it changes direction, Doppler shifts, and becomes partially po- fact that m = 1 oscillation modes are the only possible global larized. Similarly, whenever a photon is absorbed, it is not de- modesinnearlyKepleriandiscssuchasBediscs(Kato,1983). stroyed;itisreemittedlocallywith a newfrequencyanddirec- It was first applied to Be stars by Okazaki (1991) and then re- tion determinedby the local emissivity of the gas. HDUST in- visedbyPapaloizouetal.(1992). cludesbothcontinuumprocessesandspectrallinesintheopac- ityandemissivityofthegas.Sincephotonsareneverdestroyed Although the observed characteristics of the V/R varia- (absorptionisalwaysfollowedbyreemissionofanequalenergy tionswereexplainedqualitativelywellbytheglobaloscillation photonpacket),this procedureautomaticallyenforcesradiative model, it was not fully satisfactory for several reasons. First, equilibriumandconservesfluxexactly. the oscillation period is sensitive to the subtle details of disc The interaction (absorption) of the photons with the gas structureaswellasstellarparameters(Fiˇrt&Harmanec,2006). provides a direct sampling of all the radiative rates, as well Given that there are usually several stellar and disc parameters as the heating of the free electrons. Consequently, an iterative thatareonlylooselyconstrained,thishighsensitivitymeansthat schemeisadoptedinwhichtherateequationsaresolvedatthe itisalwayspossibletoobtainthesameperiodbydifferentcom- end of each iteration to update the level populations and elec- binations of parameters. For example, as far as the period is tron temperature. The process proceeds until convergence of concerned, the same effect is produced by decreasing the disc all state variables is reached. The interested reader is referred temperatureorstellarradius,orbyincreasingstellarmassorro- to Carciofi&Bjorkman (2006) for details of the Monte Carlo tationvelocity,thusmakingmodelpredictionsambiguous.Since NLTEsolution. weuseavarietyofobservationstoconstrainourmodelsinthis investigation,weareabletosignificantlyreducetheseambigui- 2.1.Geometry ties.Second,althoughthemodelpredictsoscillationmodesthat precess in the direction of disc rotation (prograde modes), as ThegeometryoftheproblemisdefinedinFig.1.Thestarrotates suggestedbyobservations(e.g.,Teltingetal.,1994),themodes counterclockwise with the rotation axis parallel to the z direc- arelessconfinedtotheinnerpartofthediscthanexpectedfrom tion,andislocatedattheoriginofthecartesiansystem(x,y,z). the observations. This can be a problem for isolated Be stars, TheCSdiscisassumedtolieintheequatorialplane(xyplane). forwhichthediscextendstoalargedistance,butweexpectthat Since we are investigating non-axisymmetric precessing discs, thisisnotsoseriousfortruncatedsmalldiscssuchasthoseinbi- we define the x-axis so thatit precesses togetherwith the disc. naryBestarslikeζTau.Asamatteroffact,wewilldemonstrate Therefore,boththexandydirectionsaretime-dependent. inlatersectionsthatalessconfinedmodebetterreproducesthe Wedefinethe(x,y,z)systemsothattheobserverislocated ′ ′ ′ observeddata. at x = and the plane of the sky is parallelto the yz plane. ′ ′ ′ ∞ This system is obtained by rotating the xyz system first by φ Inorderto criticallytestwhethertheglobaldiscoscillation degreesaroundthezaxisandthenby90 idegreesaroundthe modelandtheviscousKeplerianscenarioareapplicableforthe − y axis.Theanglei(0 i < 180 )istheviewingangleandthe CS disc of ζ Tau, it is necessary to rigorously solve the cou- ′ ≤ ◦ angleφ(0 φ<360 )describesthetime-dependentpositionof pledproblemsoftheradiativetransfer,radiativeequilibriumand ≤ ◦ thex-axis. statistical equilibriumin non-localthermodynamicequilibrium To complete the geometrical description of the system we (NLTE)conditionsforthecomplex3Dgeometriespredictedby specifytheangleγthattheprojectionoftherotationaxisonthe theglobaloscillationscenario.A3Dapproachfortheradiative plane of the sky, z, makes with respect to celestial north. We transfer is importantto account for geometricaleffects such as ′ adopt the usual convention that γ is measured east from north escape of radiation towards the optically thin polar regions, or andthatγliesintheinterval 180< γ 180 .Theanglesγ,φ theshieldingofpartofthediscbyalocaldensityenhancement. − ≤ ◦ andiareallfreeparameters. We use,forthispurpose,thecomputercodeHDUSTdescribed in Carciofi&Bjorkman (2008, 2006) and Carciofietal. (2008, 2004) 2.2.TheCentralStar InSect.2weexplainourmainmodelassumptions.InSect.3, TheadoptedstellarparametersarelistedinTable1.Asdiscussed we outlinethe adoptedcomputationalprocedureanddiscussin byRiviniusetal.(2006),determiningthepropertiesofthecen- detailsthefittingoftheobservationsreportedinPaperI.Finally, tral star is difficult because of the large Vsini (320kms 1, − inSect. 4and5we discussourresultsandpresentasummary Yangetal.,1990)andthepresenceofanopticallythickCSdisc. ofourconclusions. The adopted parameters for the central star, together with the Carciofietal.:Testingthe2DGlobalDiscOscillationModel 3 Table1.ModelParametersforζ Tau Parameter Value Type Ref.a StellarParameters SpectralType B1IVeshell fixed 2,3 R 5.9R fixed 2,4 p R 7.7R⊙ fixedb 1 e T 2500⊙0K fixed 2 p T 17950K fixedc 1 e V/V 0.8 fixed 5 crit V 530kms 1 free 1 crit − M 11.3M fixedd 1 L 7400L⊙ fixed 1 T 19370⊙K fixede 1 eff DiscParameters Σ 2.1gcm 2 free 1 0 − ρ 5.9 10 11gcm 3 fixedf 1 0 × − − R 130R fixed 1 t H 0.208⊙R fixedg 1 0 V/RP⊙arameters Fig.1.Geometryoftheproblem.Thestarislocatedatthecenter P 1429d Fixed 2 ofthe xyzsystem,withtherotationaxisalignedwiththezaxis. T0 50414h Fixed 2 Thediscisinthe xyplane.Theobserverliesalongthe x′ direc- kW2eaklineforce 50..70406 10 2(̟/R )0.1 FFrreeee 11 tion,whichisdefinedbytheviewingangleiandtheazimuthal M˙ 10 11×M −yr 1 e Free 1 angleφ,whichistheanglebetweenthexaxisandtheprojection − − α 0.4 ⊙ Free 1 of x ontothe xyplane.Theplaneoftheskycorrespondstothe ′ δ 0.95 Free 1 yz plane,andz makesanangleγwithrespecttocelestialnorth ′ ′ ′ GeometricalParameters i 95 free 1 ◦ γ 32 free 1 ◦ φ 280 free 1 0 ◦ distance 126pc fixedi 6 BinaryParameters P 132.97d fixed 7 orb e 0 fixed 7 disc parameters(Sect. 2.3), forma consistentsetofparameters a 257R fixed 7 inthesensethat,asweshallseebelow,theycanreproducewell q 0.121⊙ fixed 7 alltheobservationalproperties.However,itmustbeemphasized a References: 1) this work; 2) Paper I; 3) Yangetal. (1990); 4) that many stellar parameters are somewhat uncertain. For in- Harmanec (1988); 5) Cranmer (2005); 6) Perrymanetal. (1997); 7) stance, equally good fits are obtained by varying the physical Harmanec(1984) sizeofthestarorthephotospherictemperaturesbyabout10%. b R wascalculatedfromR andV/V forarigidlyrotatingstarin e p crit theRocheapproximation. Given the large value of Vsini, we adopt a rotationally c Theequatorialtemperature,T ,wascalculatedfromV/V andT deformed and gravity darkened central star instead of a sim- e crit p forarigidlyrotatingstarintheRocheapproximation. ple spherical star. Assuming uniform rotation, the stellar rota- d MwascalculatedfromV . tion rate is given by a single parameter V/V , i.e., the ratio crit crit e T correspondstotheaverageeffectivetemperatureofarotation- between the equatorial velocity and the critical speed V = eff crit allydeformed,gravitydarkenedstar. [2GM/(3R )]1/2, where R is the polar radius. We assume p p f ρ isgivenbyΣ (2π) 1/2H 1(Eq.10). the stellar shape to be a spheroid, which is a reasonable ap- 0 0 − 0− g H was calculated from Eq. 8, assuming a mean molecu- proximation for the shape of a rigidly rotating, sub-critical 0 lar weight of 0.6 and an electron temperature of 18000K after star (Fre´matetal., 2005). Given the stellar rotation rate and Carciofi&Bjorkman(2006). the polar radius, the equatorial radius, Re, is determined from h ForT wehavechosenthepeakvalueofcycleI(seePaperI). 0 the Roche approximation for the stellar surface equipotentials i d was considered a free parameter within the 113 — 148 pc (Fre´matetal.,2005).Forthegravitydarkeningweusethestan- range,correspondingtotheaccuracyinthedistancedeterminationthe dard von Zeipel flux distribution (vonZeipel, 1924), according Hipparcossatellite. to which F(θ) g (θ) T4 (θ), where g (θ) and T (θ) are ∝ eff ∝ eff eff eff the effective gravity and temperature at stellar latitude θ. The stellarrotationratewasassumedtobeV/V =0.8.Forthepo- 2.3.DiscStructure crit lar radius we have adopted the value R = 5.9R , typical for p We assume a steady-stateviscousdecretiondisc, forwhichthe anevolvedmainsequencestarofspectraltypeB1I⊙V(Harmanec surfacedensityisgivenby(e.g.,Bjorkman&Carciofi,2005) 1988;PaperI). M˙ GM 1/2 fittinTgheinvaolurdeeorftVocrbiteisntTraebplreod1ucweasthdeeteemrmisinsieodnfrloinmetphreodfialetas Σ(̟)= 3παc2s (cid:18) ̟3 (cid:19) h(R0/̟)1/2−1i , (1) (Sect. 3.2). The value of the stellar mass was then calculated where from V and R . Both M and V we obtain agree with the crit p crit valuestabulatedbyHarmanec(1988). ̟=(x2+y2)1/2 (2) 4 Carciofietal.:Testingthe2DGlobalDiscOscillationModel isthedistancefromthecenterofthestar,αistheviscositypa- Wealsoassumethatthediscisinverticalhydrostaticequilib- rameterofShakura&Sunyaev(1973),M˙ isthestellardecretion rium.Inthiscase,itcanbeshownthatthedensityofisothermal rate and R is an arbitrary integrationconstant. The isothermal discsisgivenby(Bjorkman&Carciofi,2005) 0 soundspeedisc =(k T)1/2(µm ) 1/2,wherek istheBoltzman s b H − b z2 constant,T isthegaskinetictemperature,µisthemeanmolec- ρ(̟,z) = ρ (̟)exp , (7) ularweight,andm isthemassofthehydrogenatom. 0 −2H2! H TheintegrationconstantR0isafreeparameterrelatedtothe whereρ0(̟)isthediscdensityatthemidplane(z = 0),andthe physical size of the disc. For time-dependent models, such as discscaleheightisgivenby thoseofOkazaki(2007),R growswithtimeandthusR isre- 0 0 latedwiththeageofthedisc. ̟ 3/2 H(̟)= H , (8) ζ Tau is a well-knownsingle line binary(Harmanec, 1984) 0 R ! e anditisusuallyassumedthatthetidalforcesfromthesecondary physicallytruncatethediscata radiusRt, correspondingtothe whereH0 ≡csVc−r1itRe tidal radius of the system (Whitehurst&King, 1991). We as- Wecanfindtheradialdependenceofρ0 byconsideringthat sumethatthediscissufficientlyoldthatR R ,inwhichcase Σistheintegralofρinthezdirection t 0 ≪ Eq.1canbewritteninasimpleparameterizedform Σ(̟)= ∞ρ(̟,z)dz= √2πHρ (̟) . (9) 0 R 2 Z Σ(̟)=Σ R2̟ 3/2 (R /̟)1/2 1 Σ e , (3) −∞ 0 e − h 0 − i≃ 0(cid:18)̟(cid:19) FromEqs.3to9wehave Σ R 3.5 z2 whereΣ0isaconstantthatmeasuresthedensityscaleatthebase ρ(̟,z)= 0 e exp . (10) ofthedisc,andiswrittenas √2πH0 (cid:18)̟(cid:19) −2H2! Beforedescribingtheglobaloscillationmodel,afewwords M˙ GM 1/2 Σ = (R /R )1/2 1 . (4) about the HDUST solution of the gas state variables are war- 0 3παc2 R3 ! 0 e − ranted.FromEqs.3and7weseethatthesolutionsforboththe s e h i verticalhydrostaticequilibriumandtheradialdiffusiondepends In Eq. 1 the disc densityscale is controlledbothby R and 0 onthedisctemperature.Inarecentwork,Carciofi&Bjorkman M˙. Fora truncateddisc system R isunknown,becausethe in- 0 (2008) studied the hydrodynamics of non-isothermal viscous formationaboutthediscageisdestroyedbythetidaldisruption Keplerian discs and found that the density structure of such oftheouterdiscbythesecondary.Therefore,insuchasystemM˙ discs can strongly deviate from the simple isothermal solu- cannotbe determinedobservationally,unlesssome other infor- tion of Eqs. 1 and 7. The dissimilarities between the two so- mationabouttheoutflowisavailable.Fromthecontinuityequa- lutions are most marked for the inner disc (̟ < 10R ), from e tion whence the polarization continuum and most o∼f the line flux comes. Indeed, a comparison between the emergent spectrum M˙ =2π̟Σ(̟)v (5) ̟ fortheself-consistentsolutionwithisothermalmodelsrevealed important differences between the two (see, e.g., Fig. 6 of weseethatthemassdecretionrateisrelatedbothtothesurface Carciofi&Bjorkman,2008). density andthe radialvelocity,v . Thus,in orderto determine ̟ Themodelspresentedin thisworkhavean isothermalden- thedecretionratetheobservationsmustprovideameasureofthe sity structure but are not isothermal, because for the assumed outflowvelocity. densitydistributionHDUSTcalculatesthefullNLTEandradia- This latter quantity is difficult to obtain directly from the tiveequilibriumproblem,therebyprovidingthegastemperature observations, because the outflow velocities in the inner disc asafunctionofpositioninthedisc.Carciofi&Bjorkman(2008) areprobablyseveralordersofmagnitudelowerthantheorbital demonstratedthatsuchmixedmodelsareamuchbetterapproach speeds(Hanuschik,1994,2000;Waters&Marlborough,1994). totheproblemthanapurelyisothermalmodel.Thereasonwhy However, the outflow velocity can, in principle, be determined anisothermaldensitystructurewasassumedinthisworkliesin indirectly.Asdiscussedbelow,theeigenvaluesoftheglobalos- the fact that at the moment only an isothermalsolution for the cillationmodedependontheviscosityparameterα.Thisparam- global disc oscillations is available. Lifting this inconsistency eter,inturn,setstheoutflowvelocity(thiscanbeseenbysubsti- betweenthecalculatedtemperaturedistributionandtheassumed tutingEq.3inEq.5).Therefore,ifwecansufficientlyconstrain densitywillbeleftforfuturework. the global oscillation eigenvalues (see Sect. 4), the value of α TomodeltheV/RvariationsandtheVLTI/AMBERobserva- and,asaconsequence,thedecretionratecouldbeobtained. tionsreportedinPaperI,Eq.3mustbemodifiedaccordingtothe From the binary parameters of Table 1, we find that the globaloscillationmodelofOkazaki(1997)andPapaloizouetal. Roche radius of the system is R = 144R . Here we have Roche (1992).Incalculatingthegravitationalpotential,wetakeintoac- usedtheapproximateformulafortheRochera⊙diusbyEggleton countthequadrupolecontributionduetotherotationaldeforma- (1983),whichisgivenby tionoftherapidlyrotatingBeprimary(Papaloizouetal.,1992). Wealsotakeintoaccounttheazimuthallyaveragedtidalpoten- 0.49q 2/3 R =a − , (6) tial(Hirose&Osaki,1993).Thepotentialisthengivenby Roche 0.69q 2/3+ln(1+q 1/3) − − ψ GM 1+k Ω⋆ 2 ̟ −2+q̟ 1+ 1 ̟ 2 , (11) wbihnearrey amaisssthraetiose.mSii-nmcaejotrheaxtiisdaalnrdadqius=is Map2p/rMo1ximisattehlye ≃ − ̟ 2 Ωcrit! Re! a " 4(cid:18)a(cid:19) # 0.9RRoche in circular binaries with small mass ratio (q < 0.3, where Ω⋆ is the angular rotation speed of the star,Ωcrit = Whitehurst&King,1991),weassumethatthediscaround∼ζTau 2(3R ) 1V , and k the apsidal motion constant. In this equa- p − crit 2 istruncatedatR =130R . tionthefirsttermis thepoint-masspotentialofthe Be star,the t ⊙ Carciofietal.:Testingthe2DGlobalDiscOscillationModel 5 secondtermisthequadrupolecontribution,andthethirdtermis (ω Ω)2 κ2 + v2κ2/c2 > 0 (e.g., Okazaki, 2000), which is − − ̟ s theazimuthallyaveragedtidalpotential. explicitlywrittenas The radial distribution of the rotational angular velocity Ω(̟) is derived from the equation of motion in the radial di- ω < Ω κ 1 v̟ 2 1/2 rection,andiswrittenexplicitlyas − − c ! Ω(̟) ≃ (cid:18)G̟dM3ln(cid:19)1Σ/2H1+2k2 ΩΩc̟⋆rit!2ǫ 1R/̟2e!−2− 2q(cid:18)̟a(cid:19)3 ≃ (cid:18)+G̟kM3 (cid:19)Ω1/⋆2"−212sR eddl2lnn+̟Σ1+ηǫdd2l̟nln̟Σǫ2+!(cid:18)1̟H(cid:19)v2̟+234q. (cid:18)̟a(cid:19)3 (19) +dln̟(cid:18)̟(cid:19) −η Re! # , (12) Here, the2 fiΩrsct,rits!ec(cid:18)o̟nd(cid:19), and2thir dRtee!rms a2re cdsu!eto the pressure undertheapproximationz2/̟2 1,wherewehaveincludeda gradientforce,tidalforceandthequadrupolecontributionofthe ≪ hypotheticalradiativeforceintheformof potentialaroundthedeformedstar,respectively.Thefourthterm is due to the radiation force and the last term is the contribu- GM ̟ ǫ tion due to advection. Near the star, the quadrupole term and F = η , (13) rad ̟2 R ! theradiationtermplaythemajorrole.Thus,theeigenfrequency e is mainly determined by these terms and the pressure gradient whereηandǫ areparameterscharacterizingtheforceduetoan term.Thelargerthecontributionofthequadrupoletermandthe ensemble of optically thin lines (Chen&Marlborough, 1994). radiation term and/orthe smaller the sound speed (or tempera- Thentheassociatedlocalepicyclicfrequencyκ(̟)iswrittenas ture),thehighertheeigenfrequency.Ontheotherhand,theother terms, i.e.,the tidaltermand theadvectiveterm, aswellas the dΩ 1/2 pressuregradienttermaffecttheoscillationbehaviorintheouter κ(̟) = 2Ω 2Ω+̟ " d̟!# part of the disc. The larger the contributionof these terms, the lesstheeigenmodeisconfined,andviceversa. = GM 1/2 1 k Ω⋆ 2 ̟ −2 2q ̟ 3 Next, we consider the effects of viscosity. When the mass- (cid:18) ̟3 (cid:19) − 2 Ωcrit! Re! − (cid:18)a(cid:19) ddeiucsreΣtionarreatfiexeM˙d,atnhde trhaediasulrvfaecloecidteynvsityisatprthoepoirntnioenradlitsoctrhae- dlnΣ d2lnΣ H 2 ̟ ǫ 1/2 viscosi0ty parameter α. Thus, one effect̟of viscosity is to make + 2 + η(1+ǫ) .(14) " dln̟ dln̟2#(cid:18)̟(cid:19) − Re! ) the oscillation mode less confined, in terms of the advective term. Katoetal. (1978) showed that m = 1 modes in nearly The perturbed surface density is obtained by superposing Kepleriandiscs, which are neutralif the viscosity is neglected, a linear m = 1 perturbation on the above unperturbed state becomeoverstablewhenviscouseffectsaretakenintoaccount. (Eq. 3) in the form of normalmode of frequencyω that varies Thegrowthrateisproportionaltoα(seealsoNegueruelaetal., asexp[i(ωt φ)].Forsimplicity,theperturbationistakentobe 2001). − isothermal.Thelinearizedperturbedequationsarethenobtained It is important to note that our model provides prograde asfollows, modes for a plausible range of parameters. The fundamental mode, in general, has the eigenfrequencyof the order Ω(R ) e "i(ω−Ω)+v̟dd̟#ΣΣ′ + ̟1Σdd̟ ̟Σv′̟ − i̟v′φ =0, (15) eκv(Rene)i.fFnroomradEiaqt.io(1n9f)oirtceisisobinvcilouudsetdh.aTthΩu(sR,teh)e−reκi(sRne)oiwsapyostoitirve−e- (cid:0) (cid:1) producetheobservedV/Rperiodofζ Taubyretrogrademodes, d Σ dv d c2 ′ + i(ω Ω)+ ̟ +v v 2Ωv =0,(16) whichhavenegativeeigenfrequencies. sd̟ Σ! " − d̟ ̟d̟# ′̟− ′φ It should also be noted that linear models like ours cannot i d Σ κ2 predicttheamplitudeoftheoscillations.Thus,forthepurposeof c2 +α ′ + v s −̟ d̟! Σ 2Ω ′̟ comparisonwith variousobservations,we will assume that the nonlinear perturbation patterns are similar to the linear eigen- v d + i(ω Ω)+ ̟ +v v =0, (17) functionsobtainedfromtheaboveequations,andtakethemax- " − ̟ ̟d̟# ′φ imumvalueoftheperturbedpartofthesurfacedensity,Σ/Σ,to ′ beafreeparameter,δ. where Σ is the Eulerian surface-density perturbation, and ′ (v ,v ) is the verticallyaveragedvelocityfield associated with ′̟ ′φ the perturbation. We solve Eqs. (15)-(17) with a rigid wall 3. Results boundarycondition,(v ,v ) = 0, at the inner edge of the disc ′̟ ′φ 3.1.ModelingProcedure andafreeboundarycondition,∆p = 0,attheouterdiscradius, where∆pistheLagrangianperturbationofpressure. The global oscillation model described above is not a trivial Itisinstructivetoconsidertheeffectsofsomeparameterson problemfor a radiativetransfer solver. A verylarge numberof them=1oscillationsbeforesolvingtheperturbationequations. grid cells (typically about100000)is needed to accurately de- For simplicity, we first neglect viscous terms. Then, the local scribethe̟,zandφdependenceofthedensity,ofthegasstate dispersionrelationisobtainedfromEqs.(15)–(17)as variablesandoftheradiationfield(fordetailsabouttheadopted cellstructureinHDUST,suchasradialandlatitudespacing,see (ω Ω kv )2 κ2 =c2k2, (18) Carciofi&Bjorkman, 2006). Solving the radiative equilibrium − − ̟ − s andNLTEstatisticalequilibriumproblemsforeachcellinvolves wherekistheradialwavenumber.Inthiscase,thepropagation a long series of iterations to accurately determinethe coupling region, where k is a real number, is given as the region where between radiationand state variablesin the dense environment 6 Carciofietal.:Testingthe2DGlobalDiscOscillationModel oftheinnerdisc.This,andthesubsequentcalculationoftheob- servables, takes between 2 and 3 days for a single model in a clusterof50PentiumIVcomputers. This is prohibitivelylong for the iterative procedureneces- sary to fit the observations. For this reason, we have initially modeledtheglobaldiscpropertiesusingasimple2Dviscousde- cretionmodelusingEqs.(3)and(7).This2Danalysisallowed us to place initial constraints on several parametersof the sys- tem,includingΣ ,V ,i,andd(Sect.3.2). 0 crit Once a suitable 2D modelwas found,we usedit asa start- ing point for the detailed 3D modeling described in Sect. 2.3. With this 3D model we performed a simultaneous analysis of theVLTI/AMBERinterferometryandtheV/Rpropertiesofthe HIlines,thatallowedustoconstrainwelltheglobaldiscoscil- lationsparameters(Sect.3.3). 3.2.Two-DimensionalViscousDecretionDiscModel As discussed in Paper I, ζ Tau has shown little or no secular evolution since the onset of the current V/R activity started in 1992, and it is argued that the global disc properties, basically itsdensityscale,havebeenapproximatelyconstantintheperiod. Therefore,modelingthedatawitha2Dviscousdecretiondiscis perseausefulexercisethatcanprovideinsightintotheaverage propertiesofthisdisc. Initially, a set of observationsthat are representativeof the period from 1992 through 2008 must be defined. For both the spectral energy distribution (SED) and polarization in the 0.3 —1µm region,we used the PBO spectropolarimetricobserva- tions made in March, 1996 (Woodetal., 1997). The flux and polarizationlevelswerescaledinordertomatchtheaverageV- band flux and polarization of the period ( V = 3.02 0.07, h i ± p = 1.46 0.08%, see Paper I). To extend the SED to the V h i ± IR, we used archival data from the IRAS point source cata- log(Beichmanetal.,1988)andthe2MASScatalog(Cutrietal., 2003). Becausethelineprofilesarevariableandstronglydependent on the disc asymmetries caused by the global oscillations, we have fitted only the average Hα peak separation in this period (240kms1,PaperI)andtheaverageHαEW( 15.5 0.1Å,Paper − ± I). Fig.2. Emergent spectrum of ζ Tau. The dark grey lines and One of the mostbeautifulcharacteristicsof the viscousde- symbolsaretheobservationsandtheblacklinesrepresentthe2D cretiondiscmodelisitsrelativelysmallnumberofparameters. modelresults. Top: Visible SED(datafromWoodetal., 1997, Todescribethesteady-statestructureofthediscweneedtospec- scaledinfluxtomatchtheaverageV-bandmagnitudefrom1992 ifythedecretionrate,M˙,theviscosityparameter,α,andtheage to the present).Middle: IR SED (data from2MASS – squares, ofthediscthatisenclosedintheparameterR (Eq.1).Wefur- 0 andIRAS–circles).Bottom: Continuumpolarization(datafrom therneedtospecifythestellarcriticalvelocity,V ,whichisa crit Woodetal., 1997, also scaled in level to match the average V- measureofthestellarmassandaffectsthediscstructurebyset- band polarization from 1992 to the present). In the two upper tingitsscaleheight(Eq.8)androtationspeeds.Inthecaseofthe panels,thelightgreylinecorrespondstotheunattenuatedstellar discaroundζ Tauweassume,asexplainedinSect.2.3,thatthe SED discistruncatedbyacompanionstar.Thismechanismerasesthe informationaboutthediscageandmakestheproblemundeter- mined.ThisbasicallymeansthatM˙,αandR areallcontainedin 0 asingleparameter,Σ ,whichdeterminesthediscdensityscale. by the Hipparcos satellite (Perrymanetal., 1997), models that 0 Thereforethediscisdescribedbyonlytwofreeparameters,Σ didnotproducedinthisrangewerediscarded. 0 andV . In additionto thesetwo physicalparameterswe must Figs. 2 shows our best fit to the data using the 2D viscous crit specifytwogeometricalparameters,theviewingangle,i,andthe decretion model. This best fit was achieved for the following distancetothestar,d. parameters For this analysis severaltens of models were run, covering a large range of values for Σ , V and i. For each model, the Σ =1.7gcm 2, 0 crit 0 − valueofdwascomputedinordertomatchtheobservedfluxlev- i=85 , ◦ els.Sincedisknowntobeintherangebetween113—148pc, V =530kms 1, crit − whichcorrespondstotheaccuracyinthedistancedetermination d=126pc, Carciofietal.:Testingthe2DGlobalDiscOscillationModel 7 a little overestimated or the geometry of inner disc is not well describedinourmodels(seeSect.4). The predicted polarization matches the observations well. The size of the Balmer jump and the level and slope of the polarization in the Paschen continuum show good agreement. However,intheBrackettcontinuumourmodelsoverestimatethe polarizationbyapproximately10%. We obtain an inclination angle of i = 85 for our best-fit ◦ model. The fact that ζ Tau is a well-known shell star (Porter, 1996; Riviniusetal., 2006) supportsour findings. In Fig. 3 we compare the model polarization for 3 different inclination an- gles, 83 , 85 and 87 . Differences in i as small as 2 produce ◦ ◦ ◦ ◦ significantchangesin the polarization,which indicatesthatthe inclinationangleiswellconstrainedbyouranalysis. TheextentthattheCSdiscmodifiesthestellarradiationcan be assessed by comparing the model SED with the unattenu- Fig.3. Continuumpolarizationforζ Tau.Themodelresultsfor atedstellarradiationinFig.2.Sincethestarisviewedcloseto three differentviewing angles, as indicated,are comparedwith edgeon,thiscausesextinctionofUV andvisibleradiationthat theobservations isreemittedatlongerwavelengths.Mostofthereprocessedra- diationescapesinthepolardirection,duetothefactthatthever- ticalopticaldepthismuchsmallerthantheradialopticaldepth. in additionto the fixed stellar parameterslisted in the first part However,aportionofthisradiationwillescaperadiallyandwill ofTable1. We recallthatforallmodelsthedisc wastruncated contributetothelargeIRexcessobserved( 3magat25µm). ≈ at130R . The Hα peak separation and EW for the best-fit model are 10−1T1hgec⊙mab−o3.veThveanluuembfoerrdΣen0sictyororefsppaorntidcslestocaρn0be=estim5.a9te×d 2av4e0rkamges−v1alauneds−fo1r4tÅhi,sretismpeecptievreiloyd,.inWgeoohdavaegrceheomseenntVwcriitthatshea assuming a mean molecular weight of 0.6, typical for ionized free parameterand,therefore,no a prioriassumption aboutthe gaswithsolarchemicalcomposition(recallthatourmodelshave stellar mass have been made. Vcrit is well constrained by the onlyH).Weobtainavalueof5.9 1013cm 3,whichischaracter- model, as both the CS emission lines and the continuum po- − isticofBestarswithdensediscs×,suchasδSco(Carciofietal., larization depend on this value (recall that Vcrit affects the disc 2006). scaleheight). Vcrit and the assumed Re = 7.7R give a stellar Beforediscussingourresults,itisusefultorecallthephysi- massof11.3M ,inagreementwithHarmanec(⊙1984).Thisfur- calprocessesthatcontroltheemergentspectrum.Inthevisible, thersupportsth⊙evaluesoftheparametersadoptedforζ Tau. the SED dependson the spectral shape of the stellar radiation, The good fit that was obtained using a relatively simple which,forarotatingstar,iscontrolledbythelatitude-dependent modelprovidesstrongobservationalsupporttoourassumptions. effectivetemperature.Thestellarradiationismodifiedbyseveral The power-law description for the surface density proved suc- processesintheCSdisc:electronscattering,bound-freeabsorp- cessfulindeterminingtheslope ofthe IRSED, andthehydro- tionandfree-boundemissionbyHIatoms,free-freeabsorption static equilibrium assumption successfully determined the cor- andemissionbyfreeelectrons,andlineabsorptionandemission rect shape and level of the continuum polarization. Therefore, byHandotherelements.Theabsorptionandemissionprocesses weconcludethatasteady-stateviscousdecretiondisc,truncated dependontheatomicoccupationnumbersandelectrontempera- atRt = 130R ,isagooddescriptionfortheaverageproperties ture,which,inturn,dependontheradiationfieldinanon-linear oftheCSdisc⊙aroundζ Tau. andcomplexway. IntheIRtheSEDiscontrolledmainlybythefree-freeemis- 3.3.GlobalOscillationModel sion of the free electrons. The continuum polarization is pro- ducedbyelectronscattering,butisalsomodifiedbytheabsorp- WenowpresentouranalysisoftheVLTI/AMBERobservations tionandemissionprocessesinthedisc. and the V/R properties of ζ Tau using the global oscillation Theabsorption,emissionandscatteringoftheradiationde- model. In additionto the parametersΣ , V , d, and i, that we 0 crit pends also on the geometry and velocity structure of the CS havediscussedabove,wenowconsidertheparametersk ,δ,M˙, 2 material.Forinstance,thelinearpolarizationisstronglydepen- α and the weak line force, so that we can consider the m = 1 dent on the number density of electrons, the disc flaring and density perturbation within the disc. Also, the time-dependent scaleheight. The slope of the IR SED, on the other hand, de- position of the m = 1 perturbation pattern with respect to the pendsmainlyontheradialprofileofthedensity,butalsoonthe directionoftheobserver,describedbytheparameterφ,mustbe discflaringandtemperaturedistribution(Carciofi&Bjorkman, determined(seeFig.1). 2006). DifferentaspectsofoursolutionareillustratedinFigs. 4to In view of the complex dependence of the emergent radia- 9, and the parameters of the best fit obtained are summarized tion on the stellar and disc parameters, the overallgood match in Table 1. Before continuingto describe our results, some ex- between the modeland observationsis remarkable.The model planations about our modeling procedure are warranted. For a reproducestheSED longwardoftheBalmer jump(3646Å)all given model, HDUST computes the emergent spectrum for an thewaydowntothelongestIRASwavelength.Theonlypartof observerwhosepositionisspecifiedbyiandφ.Tosimulatethe theSEDwhichisnotwellreproducedisthesizeoftheBalmer temporaldependenceoftheobservablesasthedensitywavepre- jump.Inourmodels,thesizeofthejumpistoolarge,meaning cesses aroundthe star,we compute,fora giveni,the emergent thateithertheamountofneutralhydrogeninthelineofsightis spectrumforseveralvaluesofφ.Theglobaloscillationparame- 8 Carciofietal.:Testingthe2DGlobalDiscOscillationModel with the predicted values (see points for MJD 53700). The ≈ second half of the cycle is better covered. For Brγ, the model seems to reproducethe data, but the V/R valuestend to cluster abovethepredictedcurve,withtheexceptionofthetwopoints forτ>0.8.ForBr15,allthedatapointsareabovethepredicted valuesinthesecondhalfofthecycle. One interesting feature of the observations that is clearly seen in Fig. 4, are the phase differences between the cycles of different lines. These phase differences have been already ob- servedinζ TaubyWisniewskietal.(2007),whosuggestedthat theymightbeassociatedwiththe helicityofthe globaloscilla- tion. The idea behindthis suggestion is that the more optically thin infrared lines form closer to the star in the disc, whereas the Hα line forming region extends further out into the disc. Becauseofthestronghelicityofthedensityperturbationpattern Fig.4. log(V/R) vs. time for the Hα (black), Brγ (green) and (see, e.g.,Fig. 5), the averageazimuthalmorphologyofthe in- Br15(red)lines.Thepointscorrespondtotheobservations,and nerdiscisquitedifferentthanthemorphologyoftheouterdisc. the lines to our best fitting model. The red dashed line corre- Therefore, the observed phase differences would be a result of spondstotheBr15V/Rcurveartificiallyshiftedinphaseby-0.2 thedifferentaveragemorphologyofeachlineformationsite. tomatchtheobservations ThephasedifferencesbetweentheHαandBrγlinesarewell reproducedbyourmodels,whichprovidesconvincingevidence tersusedinthesimulationarethenevaluatedbycomparingthe insupporttoaspiralmorphologyofthedensityperturbationpat- theoreticalamplitudeandshapeoftheV/Rcycletotheobserva- terninthediscofζ Tau.ThephaseoftheBr15lineisnotwell tions. reproduced,however.Thisline is formedin the veryinnerpart Thetime-dependentpositionofthedensitywaveisrelatedto of the disc and it seems, therefore, that the global oscillation theV/Rphaseasfollows.AccordingtotheconventionofPaperI, modeldoesnotdescribecorrectlytheshapeofthespiralpattern wedefinetheV/Rphaseτsuchthatτ=0attheV/Rmaximum. inthisregion.InFig.4,thedashedredlineshowsthemodelV/R The ephemeridesof the V/R cycle,givenby the period, P, and curveofthe Br15lineshiftedin phaseby-0.2.Thisartificially themodifiedJuliandateforτ=0,T0,weredeterminedinPaper displacedcurvereproducestheobservations.Thismayindicate IandaregiveninTable1.Weintroducetheparameterφ0,which thatthemodelpredictsthecorrectamplitudeoftheoscillations gives the position of an arbitrary point of the spiral pattern at intheinnerdisc,eventhoughtheshape,whichsetsthephase,is τ = 0. We chose this point to be the minimum of the density notcorrect. perturbationpatternatthebaseofthedisc.Theangleφisrelated tothephaseby As mentioned above, the global oscillation model failed to reproduce the Hα V/R cycle for 0.5 . τ . 0.8. However, this φ φ τ=1 − 0. (20) may not be a problem with the global oscillation model, since − 2π the 0.5 . τ . 0.8 phase interval is where the complex triple- (Notethattheangleφ decreaseswithtimefortheprogradeos- peakedHαprofilesbecomemoreprominent.InPaperI,wehave cillation mode assumed here). Since the phase is related to the presentedevidencethatthetriple-peakanomalyarisesfromthe time by τ = (T T )P 1, we find the followingexpressionfor outermostpartofthedisc,andthatitisprobablyaline-of-sight 0 − thetemporaldep−endenceofφ effect owing to the large inclination angle of the system. The failuretoreproducethispartoftheV/Rcycleislikelyduetothe φ=2π 1 T−T0 +φ . (21) factthataphysicalprocessismissinginouranalysis. 0 − P ! Fig. 5 illustrates the properties of our solution for the disc Thisequationgivesthepositionoftheminimumofthedensity global oscillation. The oscillation mode is subject to two very patternatthebaseofthestarasafunctionoftime.Foragiven stringent constraints: it must have a precessing period equal to model,avalueofφ0isobtainedbyfittingtheV/Rphaseandthe theobservedV/Rperiod(1429days,seePaperI)anditmustbe interferometricobservations. abletoreproducetheobservedamplitudeoftheV/Rvariations. Thislatterrequirementisassociatedwiththedegreeofconfine- mentofthemode,whichessentiallymeanshowfaroutintothe 3.3.1. ConstraintsfromtheV/RObservations disctheperturbationgoes.InFig.5wecomparethetheoretical In Fig. 4 we show our best fit to the V/R cycle of the Hα, Brγ V/R curve for two modes with the same precession period but and Br15 lines. For Hα, the model reproduces the first part of different degrees of confinement. The model shown in the top the V/R cycle (0 < τ < 0.5), including the shape of the curve left is the best fit modelof Table 1. Thismodelcorrespondsto andtheobservedV/Rmaximaandminima.Inthesecondhalfof a mode with a small degree of confinement, in which the os- thecycle,between0.5<τ<0.8,theobservedV/Rratiodeparts cillations extendall the way down to the outer rim of the disc. fromthequasi-sinusoidalshapethatischaracteristicofthefirst Theothermodelisshownasanexampleofamodethatalsohas half(PaperI),andthisisnotreproducedbythemodel. thesameprecessionperiodbutitismuchmoreconfined.From ThebehavioroftheV/RpropertiesoftheIRlinesaresimilar thecomparisonofthepredictedV/Rratiosforeachmode(bot- tothatofHα,butthesmallernumberofobservationsmakethe tompanelofFig.5)wesee thatthislattermodecanbereadily results less conclusive. Unfortunately, there are only three ob- discarded on the basis that it cannot account for the observed servationsinthefirsthalfofthecycle,buttheirV/Rratiosagree amplitudeoftheV/Rvariations. Carciofietal.:Testingthe2DGlobalDiscOscillationModel 9 Thefitsofthevisibilitiesandphaseswerecarriedoutinthe followingway.Foragivendiscmodelandinclinationangle,we computedsynthetic images for several wavelengthsaround the Brγline(seeFig.7)forseveralvaluesofφ,correspondingtodif- ferentV/Rphases.Thissetofsyntheticimageswerethenrotated forseveralvaluesofγ,correspondingtodifferentorientationsof thediscintheskyplane.Next,thevisibilitiesandphasesforeach ofthemodifiedimageswerecomputedandcomparedtotheob- servedvisibilitiesandphases.Aχ2minimizationprocedurewas carried out that producedthe values of γ and φ that best fitted thedata2.Finally,wemanuallyexploredtheparameterd inthe range113 148pc.However,wehavealwaysfoundthatthishad − onlyamarginalimpactonourfitand,therefore,wehavechosen tousethedistancedeterminedfromthe2Danalysis,d =126pc. Thisprocedurewasrepeatedforseveraldiscmodelsandin- clination angles. An analysis of the flux, polarization and V/R propertieswasalsocarriedoutforthesamemodels.Themodel parametersof Table 1 correspondto the valuesthat best repro- ducedalldataanalyzed,includingtheinterferometry. Our data analysis indicates that the interferometry is much moresensivitetosomemodelparametersthanothers.Themost sensitiveparameter,andhencethebestconstrained,isundoubt- edly the disc orientationγ = 32 5 . As discussed in Paper I, ◦ ± thisvalueagreeswellwithpreviouspolarimetricandinterfero- metricdeterminations.Inparticular,itiscompatiblewiththere- Fig.5. Comparison of the results for modes with different de- centlypublishedvalueγ=37 2 ,obtainedintheK bandwith ◦ ′ grees of confinement. The top left plot shows the density per- ± theCHARAarray(Giesetal.,2007).Additionally,asshownin turbation pattern for the V/R parameters of Table 1. The pic- the middle-rightpanelof Fig. 6, the overallsize of our models turedepictsthediscasseenfromaboveanddirectionofpreces- alsoagreeswiththatmeasuredbyGiesetal.(2007).Wenotethe sionofthedensitywaveiscounter-clockwise.Amoreconfined CHARAandVLTI/AMBERestimatesareindependentbecause mode, with a weak line force of 0.04(̟/R )0.1, α = 0.3 and eps Giesetal.(2007)usedabsolutelycalibratedbroadbandvisibil- k = 0.006,isdepictedinthetoprightpanel.Thebottompanel 2 ities, whereas we have used differential phases and visibilities showsthepredictedamplitudeofV/Rratioforeachmodel(solid acrossaline. line:lessconfinedmodel;dottedline:moreconfinedmodel). The χ2 map of Fig. 6 shows that the parameter φ is not as well constrained. From the fit we find that the position of the densitywaveatthetimeoftheobservationisφ(T =54081.3)= 3.3.2. ContraintsfromtheAMBER/VLTIObservations 90 25◦. The large uncertainty comes from the combined ef- ± fects of the high inclination of the disc and the fact that the Perhapsthe strongestevidencein supportof the globaloscilla- density perturbations are spread over a large range in azimuth tionscenariocomesfromthe verygoodfitof theinterferomet- (seeFig.7).UsingEq.21wefindthatthepositionofthedensity ricdata.Theformalfitoftheobservedvisibilitiesandphasesis waveatτ=0,asgivenbytheinterferometry,isφint =295 25 . shownin Fig.6forthreebaselines(markedasA, BandC,see Despitethislargeuncertainly,thisvalueisinagre0ementwi±thou◦r Fig. 7 for their orientation with respect to the disc). In the fig- best-fitvalueofφ =280 obtainedprimarilyfromthefitofthe 0 ◦ ure,theleftcolumnshowsthephasesandthemiddlecolumnthe V/Rphases.Weemphasizethatspectroscopicandinterferomet- square visibilities for each baseline. Because of the poor abso- ricdeterminationsoftheφ arecompletelyindependent,making 0 lutecalibrationofourdataset,bothobservableswerenormalized theaboveagreementrelevant. to the value obtained in the continuum. The difference in am- The analysisof the interferometricdata for differentvalues plitude of the differentprofilesare explainedby boththe base- of i (not shown) largely confirms the value for the inclination linelengths(93m,53mand130m,respectively)andorientation angle obtained from the 2D model (Sect. 3.2). In addition, in- (43◦, 99◦ and 63◦ East fromNorth, respectively).For instance, terferometryhasallowedustoliftahistoricaldegeneracyinthe baselineA correspondsto an orientationclose to theminimum determinationofi.RecallfromFig.1thattheinclinationangle elongationaxisofthedisc,whichexplainsthesmalldropinvis- goesfrom0to180 .Fori<90 thetopofthediscisfacingthe ◦ ◦ ibility.Thequalityofthefitisquantitativelyexpressedbytheχ2 observerandthecontraryistruefori>90 . ◦ mapshownintherightcolumn. In the case of an axisymmetric disc, it is really difficult to For illustration purposes we have converted the three in- distinguishbetween i > 90 and i < 90 . For an infinitelythin ◦ ◦ terferometric phases into an astrometric displacement for each disc this distinction is impossible, since the disc aspect would spectralbinacrosstheline(top-rightpanelofFig.6).Openand beidenticalforthetwocases.Foradiscofnon-zerogeometri- filled symbolsindicate the displacementalongthe disc at min- cal thickness,radiativetransfereffects create subtle differences imum and maximum elongationaxes, respectively.The typical in the aspect of each disc, but the differences are probablytoo S-shape astrometric signature of a rotating disc is clearly dis- torted:theredshiftedpartofthediscextendsfurtherawayfrom 2 For the χ2 calculation only the spectral channels in the range thecentralstarthantheblueshiftedpart(300µasand150µasre- [ 700,+700]kms 1 wereused. Thevalueofχ2 depends onthenum- − spectively).Thisisadirect,illustrativeproofofthepresenceof b−erofspectralchannels,andinthecontinuumthesechannelsarealways anasymmetryinthedisc. wellmatchedandthereforetendtoartificiallyreducetheχ2. 10 Carciofietal.:Testingthe2DGlobalDiscOscillationModel ] 10 1.2 as 0.4 m Φ [deg]a−1 00 2V [0−1]a 001...680 Astrometry [−000...022 0.4 −1000 0 1000 v [km/s] −1000 0 1000 −1000 0 1000 10 1.2 omas] 2 rm Φ [deg]b−1 00 2V [0−1]b 001...680 Diameter fCHARA [ 01 −50 0 50 100 150 0.4 Baseline Angle [deg] −1000 0 1000 −1000 0 1000 10 1.2 60 χ2min = 1.4 24 deg] 0 0−1] 1.0 g] 40 17 Φ [c−10 2V [c 00..68 γ [de 20 10 3 0.4 0 −1000 0 1000 −1000 0 1000 0 100 200 v [km/s] v [km/s] φ [deg] Fig.6. FittingoftheinterferometricobservationsforthespectralregionaroundtheBrγline.Inallpanelsthesolidlinerepresents our best-fit model of Table 1 and the dashed line a model with i = 85 . Left and Middle: Differential phases and differential ◦ visibilitiesobtainedwithAMBER/VLTI,forthreeinterferometricbaselines.ThebaselineorientationsareindicatedinFig.7.Top- Right:Interferometricphasesconvertedintoastrometricshift(photocenterdisplacement)vs.wavelength.Openandfilledsymbols representdisplacementalongthediscminimumandmaximumelongationaxisrespectively.Middle-Right:Angulardiameterderived fromaGaussianfitoftheCHARAK bandcontinuumobservationsofGiesetal.(2007).Bottom-Right:χ2mapshowingthefitfor ′ theγandφanglesforthei=95 model,withcontoursplottedforχ2 =4andχ2 =6 ◦ smalltobedetected.However,foranon-axisymmetricdiscthis thatthebottompartofthediscofζTaufacestheEarth.Notethat distinctionis possible,providedtwo conditionsare fulfilled:1) thisdegeneracyisnearlyimpossibletobeliftedwithbroadband spatiallyandspectrallyresolvedobservationsareavailableinor- interferometrydataonly,suchastheCHARAdatashowninthe dertodeterminetherotationdirectionofthediscintheplaneof right-middlepanelofFig.6.Webelievethatthisisthefirsttime the sky, so that the northward direction can be defined, and 2) thedegeneracyinthedeterminationofiwasliftedforaBestar. thediscasymmetryischiral,i.e.,thediscisnotsuperimposable onitsmirrorimage.Wenotethatweadopttheusualconvention 3.3.3. GeometricalConsiderations thattheangularmomentumvectorpointsnorth. Both of the above conditions are fulfilled for ζ Tau. The ThegeometricalpropertiesofζTauderivedaboveareillustrated VLTI/AMBER observationsunambiguouslydetermine that the in Fig. 7. In the top left panel we plot the density perturbation western side of the disc approachesthe Earth. Thus, according pattern,asseenfromabovethedisc.Thefiguremarkstheposi- tothedefinitionofnorthandthevalueofγderivedfromthein- tionoftheobserver’slineofsightfor5V/Rphases.Theτ=0.57 terferometry, we know that the north of the system is oriented line of sight correspondsto the time of the interferometric ob- 32 east of celestial north. Condition 2 is also fulfilled, since servations. The top center panel is a sketch of the disc density ◦ the disc is highly non-axisymmetricbecause of the spiral den- projected onto the skyplane at the time of the VLTI/AMBER sityperturbationproducedbytheglobalm=1oscillation.From observations. The synthetic image for the IR continuum (top thedetailedmodelfittingoftheinterferometricobservations,we rightpanel)mimicstheprojecteddiscdensityasexpected.The findthatthemodelwithviewinganglei=95 (χ2 =1.4)fitsthe bottompanelsofFig.7showmodelimagesfordifferentwave- ◦ observationsbetter than the correspondingmodelwith i = 85 lengthsacrossthe Brγ line. Thecombinationofthe discasym- ◦ (χ2 =2.4,seealsosignificantdifferencesbetweenthemodelsin metries and the effects of rotation generates a very complex Fig.6).Therefore,theVLTI/AMBERdataallowustodetermine brightnessdistribution.
The list of books you might like
