Astronomy&Astrophysicsmanuscriptno.ms (cid:13)c ESO2012 January24,2012 Tidal effects on the radial velocity curve of HD77581 (Vela X-1) G.Koenigsberger1 E.Moreno2 andD.M.Harrington3 1 InstitutodeCienciasF´ısicas,UniversidadNacionalAuto´nomadeMe´xico,Cuernavaca,Mor.62210,Mexico e-mail:[email protected] 2 InstitutodeAstrono´ıa,UniversidadNacionalAuto´nomadeMe´xico,04510,Mexico e-mail:[email protected] 3 InstituteforAstronomy,UniversityofHawaii,2680WoodlawnDrive,Honolulu,HI,96822 e-mail:[email protected] 2 Received;accepted 1 0 ABSTRACT 2 Context.ThemassoftheneutronstarinVelaX-1hasbeenfoundtobemoremassivethanthecanonical1.5M .Thisresultrelies n ⊙ on the assumption that the amplitude of the optical component’s measured radial velocity curve is not seriously affected by the a J interactionsinthesystem. Aims.Ouraimistoexploretheeffectontheradialvelocitycurvecausedbysurfacemotionsexcitedbytidalinteractions. 3 Methods.Weuse acalculation fromfirstprinciples that involvessolving theequations of motion of aLagrangian gridof surface 2 elements.Thevelocitiesonthevisiblesurfaceofthestarareprojectedalongtheline-of-sighttotheobservertoobtaintheDoppler shifts which are applied to the local line-profiles, which are then combined to obtain the absorption-line profile in the observer’s ] R referenceframe.Thecentroidoftheline-profilesfordifferentorbitalphasesisthenmeasuredandasimulatedRVcurveconstructed. Modelsarerunforthe“standard”(vsini=116km/s)and“slow”(56km/s)supergiantrotationvelocities. S Results.Thesurface velocity field iscomplex and includes fast, small-spatial scale structures. It leadsto strong variability inthe . h photospheric line profiles which, in turn, causes significant deviations from a Keplerian RV curve. The peak-to-peak amplitudes p ofmodel RVcurvesareinallcaseslargerthantheamplitudeof theorbital motion. KeplerianfitstoRVcurvesobtained withthe - “standard” rotation velocity imply mns ≥1.7 M⊙. However, a similar analysis of the “slow” rotational velocity models allows for o m ∼1.5M .Thus,thestellarrotationplaysanimportantroleindeterminingthecharacteristicsoftheperturbedRVcurve. ns ⊙ r Conclusions. Giventheobservational uncertainty inGPVel’sprojectedrotationvelocityand thestrong perturbationsseen inthe t s publishedandthemodelRVcurves,weareunabletoruleoutasmall(∼1.5M⊙)massfortheneutronstarcompanion. a [ Keywords.Stars:binaries:spectroscopi;Stars:neutron;Stars:rotation;Stars:individual:HD77581 1 v 1. Introduction processes (Woosley & Heger 2007; Timmes et al. 1996). Van 9 1 denHeuvel(2004)hasarguedthattheGPVelsystem provides ThehighmassX-raybinarysystemVelaX-1(4U0900-40)con- 6 direct evidence that there is at least one group of neutronstars sistsofapulsar(P =283s)inaneccentric(e=0.09,P =8.96 4 pulse orb thatareindeedbornwithsuchalargemass. d) orbit around the B0.5-supergiant star HD 77581=GP Vel. 1. Particularinterestinthissystemissparkedbythemassderived The procedurefor determiningthe neutron star mass in bi- 0 for the neutron star, M ∼1.8 M , since this is significantly narysystemssuchasGPVelisbasedontheradialvelocity(RV) 2 ns ⊙ curve which is derived from photospheric absorption lines ob- largerthanthecanonical1.5M upperlimitthatispredictedby 1 ⊙ servedinthespectrumoftheopticalcompanion.Impliedinthis thestandardneutronstarequationofstateaswellasthatwhichis : approach is the assumption that the RV variations truly repre- v predictedfromsupernovamodelsforthenewlyformedneutron sent the orbital motion of the star. In the case of GP Vel, the i stars. X validityofthisassumptionishighlyquestionablebecauseithas Themaximumpossiblemassofaneutronstardependsonthe r long been known that the photospheric absorptions of GP Vel equationofstateoftheultradensematterinitsinterior(Lattimer a undergoprominentlineprofilevariability.Specifically,thelines & Prakash 2007; Page & Reddy 2006). For nuclear densities developasymmetriesovertheorbitalcyclewhichleadtosystem- ∼3×1014 g cm−3, the equation of state is a “soft” one, and the aticdeviationsofthemeasuredRVdatapointsfromaKeplerian maximummassis∼1.5M .Forlargerdensities,theequationof ⊙ RVcurve.Itisbelievedthattheline-profilevariabilityisassoci- stateisa“stiff” one,anddependingontheadoptedequationof ated with tidalforcing,non-radialpulsationsandother interac- state,theuppermasslimitmaybeashighas2.9M (Kalogera ⊙ tioneffectsinthebinarysystem(vanParadijsetal. 1977b;van &Baym1996).Themajorityofneutronstarsinbinarysystems Kerkwijketal.1995;Barzivetal.2001,Quaintrelletal.2003). do indeed have m ≤1.5 M (Thorsett & Chakrabarty 1999; ns ⊙ Thisraises the questionof whetherthese systematic deviations Schwab et al. 2010; Lattimer & Prakash 2010). On the other may not lead to an artificially large RV curve amplitude and, hand, the existence of m ∼2 M neutron stars is now fairly ns ⊙ hence,anoverestimateoftheneutronstarmass. wellestablished(Demorestetal.2010),althoughitisnotclear The tidal effects are clearly present, as shown by the opti- whethersuchhigh-massneutronstarsarecreatedatbirthinthe cal light curve which displays ellipsoidal variationswith a full supernovaeventoraretheconsequenceofsubsequentaccretion amplitude∼0.1magindicatingthatthestarisstronglydistorted Sendoffprintrequeststo:G.Koenigsberger by the neutron star’s gravitational field (Jones & Liller 1973; 1 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Zuiderwijketal.1977;Tjemkesetal.1986).VanParadijsetal. etal. (2007)andMorenoetal.(2005,2011).Anapplicationof (1977b)exploredtheeffectsofthedeformationofthestaronthe the model to the B-type binary α Vir (Spica) may be found in RV curves. They assumed a circular orbit and synchronousro- Harringtonetal.(2009). tationandcomputedtheshapeofphotosphericabsorptionlines It is worth noting that the TIDES model computesthe sur- assuminganon-uniformtemperaturedistribution.However,the facemotionsoftheopticalcomponentfromfirstprinciplesand rotation of GP Vel is clearly slower than synchronous, and it allowingfornon-lineareffects.However,therearecurrentlytwo cannot be assumed to have a simple Roche geometry. Indeed, simplificationsinthemethod.Thefirstsimplificationisthatthe Tjemkes et al. (1986) showed that a model for the ellipsoidal equationsofmotionaresolvedonlyforathinsurfacelayer,in- variationsbasedontheassumptionofan“equilibriumtide”does stead of for the entire star, which hasthe disadvantagethat the notadequatelyreproducetheshapeoftheobservedlightcurve. perturbationsfromdeeperlayersin the star are neglected.This They suggested that non-linear effects might be important, but includestheexcitationofnon-radialpulsationmodeswhichhave fullnon-linearcalculationsfortidally interactingstars are only been shown to cause oscillations in the RV curve (Willems & nowbecomingavailable(Weinbergetal.2011). Aerts2002).Ontheotherhand,itsresponseisrepresentativeof Wehavedevelopeda2Dcalculation,theTIDES1 code,that theeffectsthatareproducedontheouterstellarlayer,whichis providesthetime-dependentshapeofthestellarsurfaceandits the one that is most strongly affected by the tidal forces (see, surfacevelocityfieldforthegeneralcaseofanellipticorbitand forexample,Dolginov&Smel’chakova1992).Thesecondsim- asynchronousrotation.Usingthederivedvelocityfield,theline- plificationisthatmotionsin thepolardirectionaresuppressed, profilevariabilityiscomputed(Moreno&Koenigsberger1999; allowingonlymotionsintheradialandazimuthaldirections.In Moreno,Koenigsberger&Toledano2005).Themethod,though addition, we neglect the effects of possible temperature varia- limitedtotheanalysisofthesurfacelayer,providesinsightinto tionsacrossthestellarsurface. the non-lineareffects that appear in binary stars that are not in The benefits of the model are: 1) we make no a priori as- anequilibriumconfiguration.Inparticular,itallowstheanalysis sumption regarding the mathematical formulation of the tidal ofhighlyeccentricandasynchronoussystems.Here,tidalflows flowstructuresincewe derivethevelocityfield fromfirstprin- arepresentinforcingregimesthatarefarfromequilibriumcon- ciples;2)themethodisnotlimitedtoslowstellarrotationrates figurations, and linear approximationsare inapplicable. In this nor to small orbital eccentricities; and 3) it is computationally paperweusethismodeltoexploretheeffectsontheabsorption inexpensive. line-profiles produced by the tidal flows on the B-supergiant’s The parameters that are needed for the calculation are: the surface,andtheresultingdeformationoftheRVcurve. orbitalperiod,P,eccentricity,e,inclination,i, andargumentof In Section of 2 of this paperwe summarize the method for periastron,ω ;thestellarmasses,m andm ,theradiusofthe per 1 2 calculatingthetidally-perturbedRVcurves;Section3describes primarystar,R ,itsequatorialrotationspeed,v ,thepolytropic 1 rot themethodforchosingthestellarandbinaryparameters;Section index3,n,andkinematicalviscosityofitsmaterial,ν.Inaddition, 4containstheresults;andSection5liststheconclusions. the code requiresthe depthof the surface layer,dR/R and the 1 numbersurfaceelementsforwhichtheequationsofmotionare tobesolved.Thelatterisspecifiedbythenumberoflongitudes 2. Calculationoftidally-perturbedRVcurves into which the equator is divided and the number of latitudes betweentheequatorandthepolarregion.Wehavedoneathor- 2.1.TheTIDEScodecalculations ough investigation of the TIDES code behavior (Harrington et Our method consists of computingthe motion of a Lagrangian al.,inpreparation)andfindthat500partitionsintheazimuthal gridofsurfaceelementsdistributedalongaseriesofco-latitudes direction at the equator and 20 latitudes is sufficient to resolve coveringthe surfaceof the star with massm as it is perturbed thesmall-scalestructure.Thisimpliesapproximately6800sur- 1 by its companion of mass m , which is assumed to be a point faceelements4 inthesemi-hemisphereaboveandincludingthe 2 source. The main stellar body below the perturbed layer is as- equator.Sincetheaxisofstellarrotationisassumedtobeparal- sumed to have uniform rotation. The equations of motion that leltothatoftheorbitalmotion,theperturbationsonthenorthern are solved for the set of surface elements include the gravita- and southern hemispheres are assumed to be symmetric. Table tional fields of m and m , the Coriolis force, the centrifugal 1givesadescriptionoftheinputparametersandliststhevalues 1 2 force,andgaspressure.Themotionsofallsurfaceelementsare forthoseparameterswhichwereheldconstantthroughoutallthe coupled through the viscous stresses included in the equations calculationsinthispaper. of motion. The surface layer is coupled to the interior body of Theline-profilecalculationisperformedaftertheinitialtran- thestaralsothroughtheviscosity.Thesimultaneoussolutionof sitoryphaseofthecalculationhasdampeddown.Inthecaseof theequationsofmotionforallsurfaceelementsyieldsvaluesof GPVel,thesteadystateisattainedafter∼20orbitalcycles,and the radialandazimuthalvelocityfields overthe stellar surface, theline-profilecalculationisperformedat30orbitalcycles. v andv ,respectively. r ϕ Thecalculatedvelocityfieldisthenprojectedalongtheline- 2.2.RVmeasurements of-sight to the observer to obtain the Doppler shifts required to producethe integratedphotosphericabsorption line profiles. Foreachgivensetofinputparameters,theoutputoftheTIDES WeassumeaGaussianshapeforthelocalprofiles2 andalimb- codeyieldstwoline-profilefiles,onecontainingtheline-profiles darkeninglaw of the form s(θ) = (1−u+ucosθ), with u=0.6 arisinginthetidallyperturbedsurfaceandthesecondcontaining as in the Milne-Eddington approximation. Full details of the the profiles arising from the unperturbed, rigidly-rotating sur- modelare givenin Moreno& Koenigsberger(1999),Toledano face. The lines are measured to obtain their centroid. The cen- 1 TidalInteractionswithDissipationofEnergythroughShear 3 The polytropic index of a radiative supergiant star is n ≥3, 2 Adiscussionontheeffectontheresultsoftheintrinsicline-profile Schwarzschild,1958,p.258. shapeisgiveninHarringtonetal.(2009);theuseofGaussiansissup- 4 Thenumber of azimuthal partitionsdecreases fromtheequatorial portedbytheresultsofLandstreetetal.(2009). latitudetothepolarregion. 2 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Table1.DescriptionofinputparametersfortheTIDEScode. 2.4.Anoteonorbitalphases Parameter Description Valuesused In order to be able to compare the theoretical models with the M (M ) massofstar1 seeTable3 observationaldata,theconventionsforsettingorbitalphaseφ=0 1 ⊙ M (M ) massofstar2 seeTable3 need to be specified. In this paper,we define φ=0 atperiastron 2 ⊙ R1(R⊙) stellarradiusatequilibrium seeTable3 passage. In most observationalinvestigationsφobs=0 is defined Porb(days) orbitalperiod 8.964368d attheorbitallatitude90◦ fromthelineofnodes.Thisdefinition e orbitaleccentricity 0.0898 issuchthatφ =0atthemidpointoftheX-rayeclipse(Bildsten obs i(deg) inclinationoforbitalplane 78.8,85.9 et al. 1997, van Kerkwijk et al. 1995, Barziv et al. 2001), and ω (deg) argumentofperiastron 332 per periastronpassageoccursatφ =0.17.Hence,φ =φ−0.17. β asynchronicityparameter seeTable3 obs obs 0 ν(R2 /day) kinematicalviscosity 0.12,0.22 ⊙ a absorption-linedepth 0.7 n polytropicindex 3.0 k line-profilebroadeningparam 30 3. ObservationallyconstrainedparametersofGP N num.segmentsalongequator 500 az Vel N numberoflatitudes 20 lat ThespectraltypeandclassoftheopticalcomponentinGPVel isgivenbyMorganetal.(1955)asB0.5Ib.Howarthetal.(1997) troid was computed through the weighted sum over all wave- assignthesamespectraltypebutamoreluminousclass,B0.5Iae. lengths,λi,asfollows: Eals.ti1m9a7t7eas,oafsmsu1m’sinmgasi=s7li8e.8in◦)thteor2a4n.g0eM21.7(RMa⊙w(lvsaentPaal.ra2d0ij1s1e)t. ⊙ VanKerkwijketal. (1995)derivea valueforthe stellar radius, λ = Σλi(Ci−Ii)3/2 (1) R1=29.9–30.2R⊙,undertheassumptionthatthesupergiantstar center Σ(C −I)3/2 fills its effective Roche lobe at periastron. Rawls et al. (2011) i i derivedR =31.82,underthesameassumption.Thesevaluesare 1 where the summation is carried out between the two limits consistent with the radii of B0.5Ia (26–38 R⊙) stars listed in withinwhichthe absorptionline lies. Inpractice,the limitsare thecatalogueofPassinetti-Fracassinietal.(2001),butnotwith thepointswheretheabsorptionmeetsthecontinuumlevel.Iiand thoseofB0.5Ib(18-26R⊙)stars.Hence,GPVel’sradiusisquite Ciarethelineandcontinuumintensities,respectively.Thisdefi- uncertainbutmostlikelyliesintherange26–32R⊙.Itseffective nitionofthecentroidisthesameasthatusedbythe’e’function temperature,Teff=25000±1000KwasdeterminedbySadakane intheIRAF5 subroutinesplot. etal.(1985)fromtheequivalentwidthofUV Fe IVandFeIII The RVs were also measured through a Gaussian fit to the lines,andacomparisonwithotherB0-B1Ibstars. Fraseretal. line profiles, using the IDL routine GAUSSFIT. Although a (2010)determinedlog(g/cms−2)=2.90andTeff=26500Kfrom Gaussian fit is a very poor approximation to the actual shape amodelatmospherefittothespectrum. of the lines, the derived RV’s were in general very similar to Theorbitalinclinationangle dependson GP Vel’sassumed those obtained with the flux-weighted centroid method. In the radius.Giventhelargerangeofpossiblevaluesfortheradius,it remainderofthispaperweusetheflux-weightedcentroids. is, hence,rather uncertain.Italso dependson the observeddu- rationoftheX-rayeclipse,θ ,which,however,istimevariable c and energydependent(Quaintrellet al. 2003).With these con- 2.3.Semi-amplitudeoftheRVcurve siderations in mind, most authors agree that i ≥78◦. A recent Theradialvelocitiesthatwereobtainedfromthemodellinepro- analysis by Rawls et al. (2011) yields as most probable values fileswereusedtoconstructtheRVcurves.Thesecurveswerefit i=85.9◦andi=78.8◦.8 with the function that characterizesthe orbitalmotion;i.e., the Table 2 summarizes the ranges in the parameters that are functiondescribingtheKeplerianorbit6: generally adopted for the GP Vel system. The most reliable of itsparametersarethosethathavebeendeterminedfromtheX- V =V +K [ecosω +cos(θ+ω )] (2) raypulsar’sdelaytimes;specifically,P ,e,theprojectedsemi- r 0 1 per per orb majoraxisoftheneutronstarorbit,a sini,andargumentofperi- X where astron,ω .Weadopttheseasthebasicknownparametersand per keepthemfixedthroughoutthispaper. 2π a sini K = 1 (3) Giventheconsiderableuncertaintiesthatareassociatedwith 1 P (1−e2)1/2 the remaining parameters, we opted to construct sets of self- consistentparametersforconductingthenumericalsimulations. isthesemi-amplitudeoftheRVcurve,a isthesemi-majoraxis 1 Themethodforderivingtheself-consistentparametersetsisde- of m ’s orbit, V is a constantverticaloffset (the “gamma” ve- 1 0 scribedinsections3.1–3.3.Asampleofself-consistentparame- locity),andθisthetrueanomaly. tersislistedinTable3.TheTIDEScodewasrunfirstfora se- ThefitwasperformedusingthePowellminimizationmethod lectionoftheseparameters.Wethenconstructedasecondblock inanIDLscript7.Thefreeparametersforthefitwerea1siniand of parameterswith small variationsof one or more of the self- V0.ThefittedvalueofK1 wasthencomputedusingeq.(3). consistentparametersandtheTIDEScodewasrunfortheseas well. 5 ImageReductionandAnalysisFacility,distributedbyNOAO 6 Binnendijk,1960,p.149–151 7 ThePOWELLalgorithmasdescribedinPressetal.1992,Section 8 ThisisthevaluegivenintheirTable4,althoughinthetextthevalue 10.5 quotedis77.8◦. 3 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Table2.Publishedstellarandorbitalparameters. Fortherangeinm thatwasanalyzed,theaboveexpression 1 leads to a range R =27.7–28.9 R . It is important to note that 1 ⊙ Parameter Value Reference R1 istheequilibriumradiusofm1 inthe absenceofm2 andas- P(days) 8.964368 Bildstenetal.1997 sumingno rotationaldeformation.These two effects lead to an e 0.0898 Bildstenetal.1997 asphericalshapeandauniquevalueofthestellarradiuscannot aXsini(ls) 113.89±0.13 Bildstenetal.1997 bedefined. θ (◦) 30–36 vanKerkwijketal.1995 e ω (◦) 332.59±0.92 Bildsteinetal.1997 per vsini(km/s) 116±6 Zuiderwijk1995 3.3.RVcurvesemi-amplitude,K1 56±10 Fraseretal.2010 The orbital velocity semi-amplitude is given by Binnendijk K (km/s) 17–29.7 vanKerkwijketal.1995 opt 21.7±1.6 Barzivetal.2001; (1960,p.151): 22.6±1.5 Quaintrelletal.2003 logg(cm/s2) 2.90±0.2 Fraseretal.2010 K = 2πa1 sini (7) T (◦K) 26500 Fraseretal.2010 1 P(1−e2)1/2 eff 25000±1000 Sadakaneetal.1985 V (km/s) 1105 Howarthetal.1997 Writinga =a/(1+ m1),andagainusingKepler’sThirdLaw, i(∞◦) 78.8 Rawlsetal.2011 1 m2 70.1-90. Quaintrelletal.2003 sini m −2/3 M1(M⊙) 232.24-.0203.6 RvaanwKlseerktwali.jk20e1t1al.1995 K1 =(2πG)1/3(1−e2)1/2m21/3P−1/3 1+ m21! (8) 23.8+2.4 Barzivetal.2001 23.1-2−17.0.9 Quaintrelletal.2003 Determinationsof K1 since 1976lie in the range∼17–28km/s R (R ) 31.82 Rawlsetal. (seeTable1ofQuaintrelletal.2003forasummary).Morestrin- 1 ⊙ 29.9-30.2 vanKerkwijketal.1995 gentlimitshavebeengivenbyQuaintrelletal.(22.6±1.5km/s) 30.4+1.6 Barzivetal.2001 andBarzivetal.(21.7±1.6km/s).Theselimitsarederivedfrom −2.1 26.8-32.1 Quaintrelletal.2003 fitsofKeplerianRVcurvestotheobservationaldata. 3.1.Stellarmasses 3.4.Equatorialrotationvelocity The massratio,q = m /m , is constrainedbythe knownvalue There are two widely different determinations of the projected 2 1 ofa siniandKepler’sThirdLaw P2 = 4π a3. Using q = equatorialrotationvelocity,v sini.Zuiderwijk(1995)obtained m /mX = a /a , andanda = a +a = aG((m11++m2q)),with a and v sini=116±6km/sbyfittingsyntheticprofilestotheobserva- 2 1 1 2 1 2 2 1 a =a thesemimajoraxesofthetwostars’orbits, tions.Thisisconsistentwiththe114km/sobtainedbyHowarth 2 X etal.(1997).Recently,however,Fraseretal.(2010)appliedthe P Gm sin3i 1/2 Fouriertransformmethod(Gray2008)tothespectrumandob- q= 1 −1 (4) tainedvsini=56km/s.Theyattributetheadditionalbroadening 2π (a sini)3! X observedinthespectrallinesto“macroturbulence”. AdoptingP=8.964368,a sini=113.89light-seconds(Bildstenet Itisimportanttonotethat,ingeneral,thestellarrotationve- x locityplaysaveryimportantroleindeterminingthebehaviorof al.1997),andwritingm insolarunits, 1 thestellarsurface.Forexample,fasterrotationleadstoalarger q=0.22505(sini)3/2(m /M )1/2−1 (5) deformationofthestar.Inaddition,theratioofstellarangularro- 1 ⊙ tationvelocity,ω,toorbitalangularvelocity,Ω,playsacritical Thus,forafixedvalueofi,eachvalueofm definesaunique role in the time dependenceandamplitudeof the tidal forcing. 1 value of m . The range in m =22.5–24.5 M leads to a range Thus, the synchronicity parameter β = ω/Ω, which describes 2 1 ⊙ m =1.427–2.04M forthecasesanalyzedfori=78.8◦. thedegreeofdeparturefromanequilibriumconfiguration,isof 2 ⊙ importanceinthecalculation.Notethatwhenβ=1,thesystemis insynchronousrotation,whichisoneoftheconditionsforequi- 3.2.PrimaryRadius,R 1 librium.Inaneccentricbinary,Ωchangeswithorbitalphaseand The value of the B-supergiant’s radius, R , is particularly rel- thus,βisafunctionoforbitalphase. 1 evant since the tidal forces scale as R−3. There are three con- In our model, the stellar rotation velocityis specified as an 1 straints on its value. The first relies on the assumption that it input parameter through β0 = ω/Ω0, where Ω0 is the orbital is close to filling its Roche Lobe, an assumption that requires angularvelocityatperiastron.Itcanbe convenientlyexpressed knowledge of its rotation velocity and the neutron star’s mass, as, m .The secondrelieson:1)the durationoftheX-rayeclipse, 2)ntsheorbitalseparation,a,and3)theorbitalinclination,i.The β =0.02vrot/kms−1(P /days)(1−e)3/2 (9) thirdconstraintonR reliesonastellaratmospheremodelfitto 0 R /R orb (1+e)1/2 1 1 ⊙ its spectrumwhich yieldslog g, fromwhich R can be derived 1 givenknowledgeofm1.We havechosenthelatterconstraintto Thetwodifferentvaluesofvsiniimplyβ0 ∼0.3or0.6,(for fix R1 values for the calculations. With the observational con- R1 ∼28–29R⊙ andi>78◦).Weshallrefertothecaseβ0 ∼0.6as straint on log g, the value of R1 follows from the choice of m1 the“standard”caseandtheβ0 ∼0.3asthe“slow”rotationcase. usingR = Gm1 1/2,wheregisthevalueofthesurfacegravity. 1 g Withlogg=(cid:16)2.90(cid:17)cms−2fromFraseretal.(2010), R =5.8396(m /M )1/2R (6) 1 1 ⊙ ⊙ 4 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Table3.Setsofself-consistentinputparameters. m1 m2 i R vs vf βs βf K 1 rot rot 0 0 Kep 23.6 1.468 78.8 28.36 57.09 118.25 0.300 0.622 17.37 24.0 1.708 78.8 28.60 57.09 118.25 0.298 0.617 19.87 24.2 1.830 78.8 28.72 57.09 118.25 0.296 0.614 21.12 22.5 1.427 85.9 27.69 56.14 116.30 0.302 0.626 17.71 22.7 1.546 85.9 27.82 56.14 116.30 0.301 0.624 19.02 Thistablelistsasampleofself-consistentinputparametersderivedas described in Section 3. Cols. 1 and 2 list the stellar masses, in M ; ⊙ col. 3 the orbital inclination; col. 4 the stellar radius, in R , cols. 5 ⊙ and6the“slow”andthe“fast”rotationvelocities,inkm/s;cols.6and 7 the corresponding asynchronicity parameters; and col. 8 the semi- amplitudeoftheKeplerianradialvelocitycurve,inkm/s. Fig.1.Top:Shapeofthestellarsurfaceintheequilibriumconfig- uration(case065);black/whitecolorcodingcorrespondstothe minimum/maximum radius; the map is centered on 180◦ lon- gitude. The asymmetry in the tidal bulges stems from the fact thattheorbitalseparationissmallcomparedtotheradiusofm . 1 Bottom: The radius as a function of longitude at the equator. Allthe40orbitalphasesforwhichitwascomputedareplotted, showingthatnosurfaceperturbationsarepresentinthiscalcula- tion. 5 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Table4.ModelinputparametersandRVcurveamplitudes. Case m m i β dR/R R ν K K 1 2 0 1 1 fit Kep 61a 23.6 1.468 78.8 0.300 0.06 28.363 0.12 18.37 17.37 61b 23.6 1.468 78.8 0.300 0.10 28.363 0.12 20.14 17.37 61c 23.6 1.468 78.8 0.622 0.06 28.363 0.12 17.70 17.37 61d 23.6 1.468 78.8 0.622 0.10 28.363 0.22 17.26 17.37 64d 24.0 1.708 78.8 0.617 0.10 28.602 0.22 20.01 19.87 65a 24.2 1.830 78.8 0.296 0.06 28.721 0.22 22.05 21.12 65b 24.2 1.830 78.8 0.296 0.10 28.721 0.22 24.37 21.12 65c 24.2 1.830 78.8 0.614 0.06 28.721 0.22 21.56 21.12 65d 24.2 1.830 78.8 0.614 0.10 28.721 0.22 21.96 21.12 65e 24.2 1.830 78.8 0.614 0.08 28.721 0.22 21.51 21.12 65f 24.2 1.830 78.8 0.614 0.06 28.721 0.22 21.34 21.12 71a 22.5 1.427 85.9 0.302 0.06 27.694 0.12 18.62 17.71 71b 22.5 1.427 85.9 0.302 0.10 27.694 0.12 20.76 17.71 71c 22.5 1.427 85.9 0.626 0.06 27.694 0.12 17.88 17.71 71d 22.5 1.427 85.9 0.626 0.10 27.694 0.12 17.99 17.71 72a 22.7 1.546 85.9 0.301 0.06 27.817 0.22 20.17 19.03 72b 22.7 1.546 85.9 0.301 0.10 27.817 0.22 22.11 19.03 72c 22.7 1.546 85.9 0.624 0.06 27.817 0.22 19.34 19.03 72d 22.7 1.546 85.9 0.624 0.10 27.817 0.22 19.38 19.03 72e 22.7 1.546 85.9 0.296 0.10 28.3 0.22 22.45 19.03 72f 22.7 1.546 85.9 0.613 0.10 28.3 0.22 19.49 19.03 30 23.5 1.44 78.8 0.301 0.06 28.3 0.22 18.29 17.2 31 23.5 1.44 78.8 0.301 0.10 28.3 0.22 20.33 22.1 32 23.5 1.44 78.8 0.622 0.06 28.3 0.22 17.32 17.2 33 23.5 1.44 78.8 0.622 0.10 28.3 0.22 17.06 17.2 34 23.5 1.44 78.8 0.300 0.06 28.8 0.22 18.05 18.05 34b 23.5 1.44 78.8 0.300 0.10 28.8 0.22 19.90 18.05 34c 23.5 1.44 78.8 0.620 0.06 28.8 0.22 17.55 18.05 34d 23.5 1.44 78.8 0.620 0.10 28.8 0.22 17.38 18.05 35 24.0 1.74 78.8 0.298 0.06 28.6 0.22 21.10 20.4 36 24.0 1.74 78.8 0.298 0.10 28.6 0.22 22.98 20.4 37 24.0 1.74 78.8 0.615 0.06 28.6 0.22 20.64 20.4 38 24.0 1.74 78.8 0.615 0.10 28.6 0.22 20.55 20.4 41 24.5 2.04 78.8 0.609 0.10 28.9 0.22 23.50 23.4 50 22.5 1.52 90.0 0.307 0.06 27.7 0.12 19.99 19.5 51 22.5 1.52 90.0 0.307 0.10 27.7 0.12 22.25 19.5 52 22.5 1.52 90.0 0.624 0.06 27.7 0.12 19.14 19.5 53 22.5 1.52 90.0 0.624 0.10 27.7 0.12 19.35 19.5 55 23.5 1.44 85.9 0.301 0.06 28.3 0.12 18.38 18.05 56 23.5 1.44 85.9 0.301 0.10 28.3 0.12 20.60 18.05 57 23.5 1.44 85.9 0.622 0.06 28.3 0.12 17.65 18.05 58 23.5 1.44 85.9 0.622 0.10 28.3 0.12 17.29 18.05 This table lists the input parameters for the models that were computed. Cases 61a–72d were computed with “self-consistent” sets of input parameters,asdescribedinSection3.m ,m aregivenininSolarmasses;R inSolarRadii;νinR2/day;andK andK inkm/s. 1 2 1 ⊙ fit Kep 6 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 4. Results Table 4 lists the inputparametersof the modelruns performed forthispaper.Thefirstblockofmodels(cases61a–72d)arerun with self-consistent sets of parameters, derived as is described intheprevioussection.Theparametersofthesecondblock(list startingwithcase72e)arenotfullyself-consistent.Columns2– 5,andcolumn7list,respectively,m ,m ,i,β andR ;column6 1 2 0 1 containsthedepthofthesurfacelayer,dR/R ,usedinthemodel. 1 Columnn8liststhevalueofthekinematicalviscosity,ν(inunits ofR2/day). ⊙ 4.1.Theequilibriumcase It is illustrative to first analyze an equilibrium case having pa- rameterssimilartothoseoftheGPVelsystem,butinwhiche=0 and β=1. We chose the case with m =24.2 M , m =1.83 M , 1 ⊙ 2 ⊙ R =28.721 R and v =118 km/s, the latter corresponding to 1 ⊙ rot the“standard”rotationspeedandanorbitalinclinationi=78.8◦. Inorderforthissystemtobeinequilibrium(i.e.,β=1),wehad tosettheorbitalperiodto12.2days.We preferedmodifyingP insteadofv becauseinordertomakeβ =1usingGPVel’sor- rot 0 bitalperiod,wewouldneedtosetv ∼160km/s,significantly rot largerthanitsactualrotationspeed.Thisleadstoastrongerde- formationof the star. The main difference of using a largeror- bitalperiodinsteadoftheactualperiodisthatthetidalforceis weaker,sothetidaldeformationissmaller. ThestartoftheTIDEScalculationischaracterizedbyatran- sitoryphaseduringwhichthestaradjustsfromitsinitiallyspher- ical shape to the new equilibrium shape. During this transitory phase,thestarundergoesinitiallylargeamplitude(∼0.03R in ⊙ thepresentcase)oscillationsthatrapidlydampout,reachingan Fig.2. Top: Shape of the stellar surface for “standard” cases amplitude <0.002 R by 20 cycles after the start of the calcu- ⊙ (model65c)atperiastron.Color codingand orientationare the lation. At this time, the star has attained its equilibrium shape, same as in Figure 1. Bottom: The crosses represent the radius whichconsistsoftwotidalbulges.Figure1showsacolor-coded at each azimuth angle along the equator. The dotted line is the Mollweiderepresentationofthestellarradiusateachsurfaceel- equilibriumshapeshownin Figure1,rescaledbya factorof3. ement. The left edge of the map correspondsto azimuth angle The strong departures from the equilibrim configuration when 0◦ (i.e., the sub-binary longitude, defined as the longitude that β =0.6areevident.Thedetailedpatternchangesasafunctionof intersects the line connecting the centers of the two stars) and 0 orbitalphase,causingtheline-profilevariations. therightedgeisat360◦.Theplotbelowthemapillustratesthe magnitude of the radius at the equatoriallatitude as a function ofazimuth.Note thatthesize ofthe two tidalbulgesisnotthe 4.2.The“standard”case sameduetothefactthattheorbitalseparationisrelativelysmall comparedto thestellar radius.Theheightofthe primarybulge Thesituationisverydifferentwhenthesystemdepartsfromthe abovethe equilibriumradiusin thiscalculationisδR=0.06R⊙. equilibrium configuration.The fact that β0 ,1 implies that the ThiscorrespondstoδR/R =0.002,whichissignificantlysmaller star rotates asynchronously, and this introduces perturbations 1 thanthedepthofthesurfacelayerusedintheTIDEScalculation, on the stellar surface. We illustrate this case with a model for dR/R =0.06. Also, note that the primary bulge points directly m =24.2M ,m =1.83M ,R =28.721R (case65c).Contrary 1 1 ⊙ 2 ⊙ 1 ⊙ towards the companion, another indication that the system is totheequilibriumcase,thestellarsurfacedoesnotattainasim- in the equilibriumconfiguration.In non-equilibriumconfigura- ple shape with two tidal bulges. The map in Figure 2 shows tions,thebulgeeither“lags”behindorisadvancedwithrespect a color-coded representation of the stellar radius at each sur- tothelineconnectingthecentersofthestars,aphenomenonthat face element, with white/black indicating maximum/minimum leadstothetidaltorquesinthesystem. extent.The map correspondsto the time of periastronpassage. Once the equilibriumconfigurationisattained,there areno The Mollweide representation and orientation are the same as significantsurfacemotions.Figure1(bottom)isactuallyaplot inFigure1.Theprimarytidalbulgecanbeseentolie between ofthe equatoriallatitudeof thestar atthe 40orbitalphasesfor 320◦ and 360◦, “lagging” behind the line connecting the cen- whichitwascomputed.Thewidthofthecurveisanindication ters of the two stars (which lies at azumuth 0◦). This is as ex- ofthestabilityovertimeoftheconfiguration.Asaconsequence, pected for a sub-synchronously rotating star. Furthermore, the nosignificantvariabilityinthephotosphericabsorption-linepro- bulge shape is not smooth, but consists of smaller-scale struc- filesisobservedinourcalculations.9 tures.Thesecanbemosteasilyvisualizedintheplotshownbe- lowthemap,whichshowstheradiusattheequatorplottedasa 9 Notethatourmodeldoesnottakeintoaccountavariabletemper- aturestructureacrosstheobservablestellardisk.Inthecaseofamore phase-dependentvariationintheline-profilesisexpectedduetogravity complete calculation in which a model stellar atmosphere is used, a darkeningeffects,seeZuiderwijketal.1977;andPalate&Rauw,2011. 7 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Fig.3. Predicted line profiles at different orbital phases for case 65c. The dotted profiles correspond to the non-perturbed, rigidly-rotatingsurface.Theasymmetriesintheperturbedline- profilesareresponsibleforthedeparturesfromaKeplerianRV curve. functionofazimuth.Asimilarconfigurationobtainsforallother orbital phases as well. The comparison between Figure 2 and the equilibriumconfiguration shown in Figure 1 shows that an “equilibriumtide”representationforGP Velisa poorapproxi- mation. Thenon-equilibriumconfigurationleadstolarge-scaleflows onthestellarsurface(seeTassoul1987andEggletonetal.1998 fora gooddescriptionofthisphenomenon).Theseare referred to as “tidalflows”. The tidal flows are motionsof localizedre- gionsofthestellarsurfacerelativetotheunderlying,(assumed) rigidly-rotatingstellarinterior.Thechangingflow patternslead tovariableline-profiles.Asampleofthetypeofvariabilitypre- dicted by the TIDES code is illustrated in Figure 3, where the perturbedlineprofilesarecomparedtotheircorrespondingnon- perturbed line profiles. Note the appearance of “bumps” and Fig.4.Radialvelocitycurvesfromthetidallyperturbedlinepro- “wiggles” in the profiles, as well as the occassional blue or files (squares) compared to the actual projected orbital motion redextendedwing.Incalculationsperformedwithsignificantly (dotted curve) and best-fit curve (dash curve). ∆RV is the dif- smaller“turbulent”speeds(forexample,10km/sinsteadofthe ference between the perturbed and non-perturbed RV values. 30km/susedhere),the “bumps”are narrower,andthe profiles Top: mns=1.830 M⊙, β0=0.614 (the ”standard” case; case65c); give the appearance of having narrow discrete absorption fea- Bottom:mns=1.468M⊙,β0=0.300(case61a). tures that generally travel from “blue” to “red” along the line profile (see Harrington et al. 2009 for an example of this type ofbehavior).InthecaseshowninFigure3thediscreteabsorp- 1.0.Despitethe largedeformationsontheperturbedRV curve, tionsaresignificantlysmoothedoutduetothelarge“turbulent” thebest-fitKeplerianisverysimilartothecurvewhichdescribes speed. the actualorbitalmotion.Thus, althoughthe peak-to-peakam- TheRVcurvederivedfromtheselineprofiles,usingtheflux- plitudeoftheperturbedRVcurveissignificantlylargerthanthe weightedmeanmethod,isplottedinFigure4(top).10.Alsoplot- actualorbitalmotion,the“excursions”oftheRVcurvearesuch ted in this figure is the RV curve corresponding to the actual thattheynearlycanceloutinthefittingalgorithm.Thisindicates orbital motion (dots), and the best-fit Keplerian curve derived thatKepleriancurvefitstothetidally-perturbedRVcurvesyield fromthe analysisthatis describedin Section2.3(dashes).The reasonableresultsaslongastheentiresetofdataisfit;i.e.,with- deviationsoftheperturbedRVcurvewithrespecttothatofthe outexcludingportionsoftheRVcurve. actual orbital motion are shown below the RV curve. They are Thegeneralcharacteristicsdescribedaboveappearinallof aslargeas±5km/sinthephaseintervalsφ=0.1–0.3andφ=0.9– the model runs listed in Table 4 for which β0 ∼0.6. The semi- amplitudesofthecorrespondingRVcurvesderivedfromthefits 10 AsimilarcurveisderivedifweplottheRV’sobtainedbyfittinga (Klsq)areplottedasfilledsymbolsinFigure5withm2=mns on Gaussiantothelineprofiles theabscissa.Also plottedinthisfigurearethevaluesof K1 for 8 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 plitudesand, thus, it is crucialto have a well-established value forthisparameterinordertoproperlyestimatethemagnitudeof theRVcurveperturbations. 4.4.Dependenceonlayerdepth Calculationsperformedwithdifferentlayerdepthsinthe“stan- dard” rotation case yield similar results. However, significant differences are present in the “slow” rotation cases where the dR/R =0.10 layer depths lead to larger RV curve amplitudes 1 than those for dR/R =0.06. Figure 6 illustrates the radial ve- 1 locitiesobtainedfromcalculationswiththesetwolayerdepths. Alsoplottedarethecorrespondingbest-fitKepleriancurvesand the curve that describes the actual orbital motion (continuous curve). The largestamplitudesoccur for β ∼0.3 and dR/R =0.10, 0 1 asillustratedinFigure7whereweplotK asafunctionofm . lsq ns Thesemodelcalculationsindicatethattheβ ∼0.3systemswith 0 m aslowas1.47M areconsistentwiththeobservations. ns ⊙ Fig.5. Predicted semi-amplitude of the GP Vel RV curves for 4.5.The“bluedip”andastructuredwind differentassumedneutronstarmassesandmodelinputparame- AprominentfeatureinourmodelRVcurvesisa“bluedip”that ters:opensymbols:β ∼0.3;filledsymbols:β ∼0.6.Allmod- 0 0 appearsatφ∼0.3,shortlyafterthemaximuminthecurve.Itisa els were run with dR/R =0.06. The dotted lines indicate the 1 persistentfeatureinallourcomputations,anditiscausedbythe observed range of semi-amplitudes from Barziv et al. (2001) asymmetricalshapeofthelineprofiles.Figure8displaystheline and Quaintrell et al. (2003). The Keplerian semi-amplitudes profiles overthe phase interval0.12–0.36for model65c, illus- (crosses)aregenerallysmallerthanthosepredictedbythemodel tratingthesourceofthisfeature.Betweenphases0.14–0.22,the calculations. Note that for the slow-rotation (β ∼0.3) models, 0 linedevelopsanexcessblueabsorptionwingwhileatthesame masses as low as m =1.55 M yield K values within the ob- ns ⊙ 1 time the red wing becomes less extended. The combination of servedrange. thesetwoeffectsmovesthecentroidofthelinetowardsshorter wavelengths.Anaturalquestioniswhatcausesthisasymmetry theKeplerianorbit(crosses).Theobservationallimitsaredrawn inthelineprofiles? with horizontal dotted lines. Only models with m =m >1.7 Amapoftheazimuthalvelocityperturbationsoverthevis- ns 2 M have values of K that lie within the observationallimits. ible surface of the B-supergiant at orbital phase φ=0.2 is pre- ⊙ lsq Hence, we conclude that for the β ∼0.6 group of models, the sented in Figure 9. The light colorindicates motionof the sur- 0 systematicshiftsinRVcausedbytidalflowsproduceonlysmall faceelementsthatisfasterthantheunderlyingrigid-bodyrota- departuresfromtheactualorbitalRVcurve,andthederivedneu- tionrate.Hence,thelargewhiteareaontheleftsideofthemap tronstarmassisnotseverelyaffectedbythetidalflows. indicates that the surface elements that lie near the limb of the starareapproachingtheobserverfasterthanthestellarrotation velocity.This leads to the more extendedblue wing in the line 4.3.The“slow”rotationcase profiles. Also evident in Figure 9 is the dark area in the map, whichcorrespondsto surfaceelementswhoseazimuthalveloc- Theslowrotationmodelsalsoshowsignificantline-profilevari- ity is slower than that of the stellar rotation. When this region ability. The deformation of the RV curves leads to systemati- reaches the right limb of the visible disk, it leads to a less ex- callylargervaluesofK .Toillustratethiscaseweuseamodel lsq withm =23.6M ,m =1.468M ,andR =28.363R (case61a). tendedredwingintheabsorptionlineprofile. 1 ⊙ 2 ⊙ 1 ⊙ Figure 4 (bottom) shows that the perturbed RV curve departs An extensive analysis of GP Vel’s observational RV curve significantly more from its Keplerian motion than the β ∼0.6 has been made by Barziv et al. (2001) and Quaintrell et al. 0 case.Thereare∼5km/sexcessesaroundφ=0.1and0.65,phases (2003),using independentdata sets. Barziv etal. (2001)report which lie close to the extrema in the actual orbital RV curve. thepresenceofa“blueexcursion”intheHδdatathatproduces Thiscausesthebest-fitRVcurve(dashes)tohaveasignificantly alocalminimumintheRVcurveatφ ∼0.37.Theyinterpretthe larger amplitude than the curve describing the orbital motion “blueexcursion”intermsofaphotoionizationwake.Figure10 (dots). isaplotoftheHδ4102ÅlineRVvalues(crosses)fromBarzivet ThesummaryofKeplerianfitsemi-amplitudes,K ,forthe al.(2001)showingthatthe“blueexcursion”initiatesatφ∼0.25. lsq β ∼0.3casesisshowninFigure5withopensymbols.Formod- Thiscoincideswiththe“bluedip”inourmodels,butitsduration 0 elsinwhichthelayerdepthisdR/R =0.06,weseethatmasses isfarlongerthanthatofthe“bluedip”. Analternativeexplana- 1 aslowasm =1.55yieldvaluesof K (marginally)consistent tion for the “blue excursion”is the presenceof enhancedmass ns lsq withtheobservations. outflow after periastronpassage and it is tempting to speculate Thus, forthe β ∼0.3groupof models,the systematicshifts that it could be triggeredby the strongtidal effects that appear inRVcausedbytidalflowsproducesignificantdeviationsfrom afterperiastronpassage. theactualorbitalRVcurve. ItisalsointerestingtonotethatKreykenbohmetal.(2008) Thebasicconclusionofthissectionisthatdifferentrotation detected flaring activity and temporary quasi-periodic oscilla- speedsleadtotidalflowswithdifferentcharacteristicsandam- tions in INTEGRAL X-ray observations. In particular, the two 9 Koenigsberger,Moreno&Harrington:RVcurveofVelaX-1 Fig.7.SameasinFigure5,butfromcalculationsperformedwith alayerdepthdR/R =0.10. 1 Fig.6. Radial velocities obtained from calculations with dif- ferent layer depths for the “standard” rotation velocity (top; cases61c,d)andthe“slow”velocity(bottom;cases61a,b).Also plotted are the correspondingbest-fit Keplerian curves and the curvethatdescribestheactualorbitalmotion(continuouscurve). Fig.8. Montage of line profiles at orbital phases φ=0.12–0.36, CrossesandthedottedcurvecorrespondtodR/R =0.06;andtri- 1 showingthatthe“dip”intheRVcurveislargelyaconsequence anglesandthe dashcurvecorrespondto dR/R =0.10.Theper- 1 of the more extended blue wing and less extended red wing. turbationsduetotidalflowsaremoresignificantforthe“slow” Dottedcurvesaretheunperturbedlineprofiles.Phaseswithre- rotationvelocitycases. specttoperiastronpassagearelisted. largestflares observedin these data occuraroundorbitalphase 0.4 (with respect to periastron) in two different orbital cycles. Thisis significantbecausethisis just∼0.03in phaselaterthan the minimum in the “blue excursion” in Barziv et al.’s data. tronstar,thelargeraccretionratesleadtoanincreaseintheX- Assuming a mean wind velocity of 1105 km/s (Howarth et al. rayemission,asdescribedbyKreykenbohmetal.(2008).These 1997),agasstreamthatoriginatesneartheB-supergiantsurface authorsalsofoundtemporaryquasi-periodicoscillationswith a wouldtravel∼3×1012 cmduringthistimeinterval,whichison timescale of ∼2 hr which they attribute to alternating high and theorderofmagnitudeofthedistancetotheneutronstar.Hence, low density regions within GP Vel’s stellar wind. A similarly thepicturewhichemergesisoneinwhichtidalflowsmaycause structuredwindhasbeensuggestedtoexistintheBe/X-raybi- instabilities that produce time-dependent ouflows leading to a narysystem2S0114+650andtopossiblybeproducedbytidally structuredwind.Whenhigh-densitywindregionsreachtheneu- inducedpulsations(Koenigsbergeretal.2006). 10