Jet formation from massive young stars: Magnetohydrodynamics versus radiation pressure
Bhargav Vaidya, Christian Fendt, Henrik Beuther, Oliver Porth
DDraft version November 14, 2018
Preprint typeset using L A TEX style emulateapj v. 11/10/09
JET FORMATION FROM MASSIVE YOUNG STARS:MAGNETOHYDRODYNAMICS VERSUS RADIATION PRESSURE
Bhargav Vaidya , Christian Fendt, Henrik Beuther, and Oliver Porth Max Planck Institute for Astronomy, K¨onigstuhl 17, D-69117 Heidelberg, Germany
Draft version November 14, 2018
ABSTRACTObservations indicate that outflows from massive young stars are more collimated during their earlyevolution compared to later stages. Our paper investigates various physical processes that impacts theoutflow dynamics, i.e. its acceleration and collimation. We perform axisymmetric MHD simulationsparticularly considering the radiation pressure exerted by the star and the disk. We have modified thePLUTO code to include radiative forces in the line-driving approximation. We launch the outflow fromthe innermost disk region ( r <
50 AU) by magneto-centrifugal acceleration. In order to disentangleMHD effects from radiative forces, we start the simulation in pure MHD, and later switch on theradiation force. We perform a parameter study considering different stellar masses (thus luminosity),magnetic flux, and line-force strength. For our reference simulation - assuming a 30M (cid:12) star, we findsubstantial de-collimation of 35% due to radiation forces. The opening angle increases from 20 ◦ to32 ◦ for stellar masses from 20M (cid:12) to 60M (cid:12) . A small change in the line-force parameter α from 0.60to 0.55 changes the opening angle by ∼ ◦ . We find that it is mainly the stellar radiation whichaffects the jet dynamics. Unless the disk extends very close to the star, its is too small to have muchimpact. Essentially, our parameter runs with different stellar mass can be understood as a proxy forthe time evolution of the star-outflow system. Thus, we have shown that when the stellar mass (thusluminosity) increases (with age), the outflows become less collimated. Subject headings: accretion, accretion disks – ISM: jets and outflows – methods: numerical – MHD –radiative transfer – stars: formation – stars: massive INTRODUCTION
Outflows and jets are integral processes of star forma-tion. They are believed to be essential for the angularmomentum evolution of the cloud core and the protostar- either directly as stellar winds or indirectly by chang-ing the structure and the evolution of the surroundingaccretion disk. Also, outflows from young stars providean important feedback mechanism to return mass andenergy into the ambient medium from which the youngstar is born.Most of the current understanding about the forma-tion and propagation of jets/outflows comes from obser-vations of low-mass stars. In this case we are fortunateto know the leading dynamical parameters such as out-flow velocity, density, and temperature (e.g. Hartigan &Morse 2007, Ray et al. 2007). However, in case of mas-sive young stars many of these parameters are poorlyknown, although large multi-wavelength studies for mas-sive star forming regions have been done over the lastdecade (Beuther et al. 2002a, Stanke et al. 2002, Zhanget al. 2005, L´opez-Sepulcre et al. 2009, L´opez-Sepulcreet al. 2010, Torrelles et al. 2011). These studies suggestthat outflows are an ubiquitous phenomenon not only forlow-mass stars, but also in massive star forming regions.The standard framework for the launching process ofjets or outflows from low-mass protostars (and mostprobably also for extragalactic jets) is the model of a disk
Email: [email protected], [email protected] Member of the
International Max Planck Research Schoolfor Astronomy and Cosmic Physics at the University of Hei-delberg (IMPRS-HD), and the
Heidelberg Graduate School ofFundamental Physics (HGSFP) wind accelerated and collimated by magneto-centrifugaland magnetohydrodynamic forces (Blandford & Payne1982, Pudritz et al. 2007). A number of numerical sim-ulations of jet formation have been performed which allconfirm this picture of self-collimated MHD jets for lowmass stars (Ouyed & Pudritz 1997, Krasnopolsky et al.1999, Fendt & ˇCemelji´c 2002, Ouyed et al. 2003, Fendt2006, Fendt 2009). Its applicability has also been demon-strated for the case of extragalactic jets (e.g.Komissarovet al. 2007, Porth & Fendt 2010).The formation of the jets around young, still forminghigh-mass stars takes place in the deeply embedded colddust and gas cocoons exhibiting large visual extinctionof the order 100 to 1000 mag (Arce et al. 2007). Recentprogress in observations of outflows from young high-mass stars comes from mm and cm wavelengths. Earlylow-spatial-resolution single-dish studies suggested thatmassive outflows may have a lower collimation degreethan those of their low-mass counterparts (Shepherd &Churchwell 1996). Also, one of the famous outflow, theOrion-KL system, exhibits a more chaotic and not colli-mated structure (e.g., Schulz et al. 1995). However, laterstudies indicated that the collimation degree of jets fromhigh-mass protostars can be as high as from low-massregions (e.g. Beuther et al. 2002a, Beuther et al. 2002b,Gibb et al. 2003, Garay et al. 2003, Brooks et al. 2003,Beuther et al. 2004, Davis et al. 2004). Typical outflowrates in high-mass systems are estimated to be a fewtimes 10 − to 10 − M (cid:12) yr − , implying accretion rates onthe same order of magnitude (e.g.,Beuther et al. 2002b).The originally puzzling result that different studiesfound different degrees of jet collimation for high-massstar-forming regions could qualitatively be resolved when a r X i v : . [ a s t r o - ph . S R ] A ug it was realized that these studies in fact targeted differ-ent evolutionary stages. While jets found in the youngeststar forming regions appear similarly collimated as theirlow-mass counterparts, jets from more evolved (massive)stars exhibit much broader outflow cones. To accountfor this effect Beuther & Shepherd (2005) proposed anevolutionary picture for high-mass outflows, where thejet formation starts with similar MHD acceleration pro-cesses compared to the low-mass models during the ear-liest evolutionary stages. However, as soon as the centralsources gain significant mass, other processes come intoplay, for example, the radiation from the central star ormore turbulence at the base of the jet. Any combinationof this and other processes was expected to lower the de-gree of collimation. However, the evolutionary sequenceproposed by Beuther & Shepherd (2005) was largely ob-servationally motivated and could only provide a quali-tative explanation. A thorough theoretical treatment toaccount for this observational result is still missing - atopic which we now address in the present paper.How the radiation field affects the formation of a jetis not obvious a priori . In order to quantify and todisentangle the physical processes involved, a detailednumerical investigation is required. Essentially, stellar(and disk) radiative forces may affect jet acceleration andcollimation directly (neglecting ionization, heating, andprobably turbulent stirring for simplicity), but also in-directly by changing the physical conditions of the jetlaunching area, thus governing the mass loading or theinitial entropy of the ejected jet material. For example,numerical MHD simulations have shown that jets withhigher (turbulent) magnetic diffusivity are expected tobe substantially less collimated (Fendt & ˇCemelji´c 2002).Our previous studies (Vaidya et al. 2009) have shownthat the inner accretion disk around massive protostarsis ionized, has a high temperature and is gravitationallyand thermally sufficiently stable in order to provide asuitable launching area for an outflow. Together withthe recent observations of magnetic fields around highmass protostars (e.g Vlemmings 2008; Vlemmings et al.2010,Girart et al. 2009) , this indeed supports the pictureof a scaled-up version of low-mass stellar jet formation.In this paper we will present a detailed investigationof how a strong radiation field impacts the structure anddynamics of a magneto-hydrodynamical driven jet. Mo-tivated by the presence of strong jets and outflows inmassive star formation, we apply the standard picture ofMHD jet formation known for astrophysical jets and putit in the physical environment of a massive young star. MODEL SETUP: JET FORMATION FROM MASSIVEYOUNG STARS
In this section we discuss the model setup applied forour numerical study concerning the formation of jets andoutflows from massive young stars. The central point isthat we are going to consider the main features of the standard model of MHD jet formation which is well es-tablished for low-mass young stars or Active GalacticNuclei also for high-mass young stars. Our model con-sists of the following essential ingredients (see Fig. 1). • A central massive young star which is rapidly evolv-ing in mass M ∗ , luminosity L ∗ , and radius R ∗ . • A surrounding accretion disk with high accre-
Figure 1.
Sketch showing our model setup of inner regions arounda massive young star. Several constituents are considered: Thedisk outflow ( dashed arrows ) is launched along the magnetic fluxsurfaces. The inner most launching point is denoted by R in , jet .This outflow is subject to radiation forces ( white arrow ) from thestar and (potentially) from an inner hot accretion disk. The stellarradius R ∗ , could be as large as ∼
100 R (cid:12) , while the stellar magneticfield structure and strength B ∗ is rather uncertain. The locationof an inner disk radius R in , disk (if existent) is also not known. tion rate, estimated to be of the order of ˙ M (cid:39) − M (cid:12) yr − . • A jet launching inner accretion disk. The extensionof this area towards the star is not known and maydepend on the existence of a strong stellar mag-netic field and stellar radiation pressure. Insteadof introducing an inner disk radius R in , disk , we willrefer to an inner jet launching radius R in , jet , whichwe presume to be between 0.1 and 1.0 AU, andto which the length scale of the simulation will benormalized l = R in , jet . • A magnetic field around the protostellar object.Magnetic fields are essential for generating colli-mated high-speed outflows. Observational indica-tion for such fields around high-mass protostars ex-ists. • A strong radiation field of the high-mass young starwhich may influence accretion and ejection pro-cesses. In case when the accretion disk reachesdown to radii close to the stellar surface (no gapas in case of low mass stars), also the high disk lu-minosity may play a role for the outflow dynamics.In the following we briefly discuss the observationaland theoretical background of these constituents and fi-nally mention the limits of our model and possible modelextensions (a more detailed discussion is provided beforethe summary).
The central massive young star
HD jets and radiation pressure 3It has been proposed that outflows from massive youngstars may follow an evolutionary sequence such that out-flows tend to be more collimated and similar to jets fromlow-mass stars in the early stages of stellar evolution,whereas at later times the outflows are less collimated(Beuther & Shepherd 2005). Since this intrinsically cor-responds to an evolutionary sequence in stellar mass, wehave investigated simulations with different central mass,ranging from 20 M (cid:12) to 60 M (cid:12) .A higher stellar mass automatically gives rise to afaster outflow (supposed the relative launching radius isthe same), just because the outflow originates deeper inthe gravitational well.
The accretion disk and the jet launching area
Jets from low mass stars are thought to be launchedfrom less than 1 AU of the inner accretion disk(e.g. Anderson et al. 2003, Ray et al. 2007). In T Tauristars, these regions could be well studied via NIR inter-ferometry (e.g. Akeson et al. 2000, 2005, Dullemond &Monnier 2010). However, to probe regions ≤
100 AUaround a young high-mass star is difficult due to thelarge extinction. We had therefore studied this regionvia semi-analytic modeling in a previous paper (Vaidyaet al. 2009) to get handle on the physical properties ofsuch disks. We have shown that here the disk may reachhigh temperatures ∼ K leading to sublimation of mostof the dust and ionizing the bulk of the material.The inner disk radius in case of low mass stars is usu-ally estimated assuming magnetic pressure balance withthe accretion ram pressure. Typical values are ∼ − (cid:63) .For young massive stars, such an estimate is not possibledue to lack of knowledge on stellar magnetic fields duringthe formation stage. In addition, radiative force from thebright luminous massive star could also influence dynam-ics of the inner disk. Although our semi-analytic mod-eling indicated that the disk could extend right down tothe central star, detailed 3D models with accurate radia-tive treatment are required to get more insight in theseclose-by regions (Vaidya et al. in prep ).We have chosen the inner launching point at a dis-tance of l = 1 AU from the star (i.e. similar to low massstars). A value of l < . l >
10 AU will result in slow outflows, barely affectedby stellar and disk radiation forces.The density ρ at the inner launching point is used for aphysical scaling. We estimate ρ from the observed massfluxes, which are typically of the order of 10 − − − M (cid:12) yr − , providing ρ (cid:39) − − − g cm − ( § The magnetic field
For many years the role of magnetic fields in massivestar formation was not really known. However, recentobservations have detected relatively strong magneticfields in massive star forming regions (Vlemmings 2008).Polarimetric observations of the hot massive molecularcore HMC G31.41 have revealed a large-scale hourglass-shaped magnetic field configuration (Girart et al. 2009). Beuther et al. (2010) detected a magnetic field alignedwith the molecular outflow via polarimetric CO emis-sion. Further observations have detected synchrotronemission from the proto-stellar jet HH 80-81, indicatinga ∼ . ∼
10 M (cid:12) is in the range of massive stars(Carrasco-Gonz´alez et al. 2010).Vlemmings et al. (2010) have measured magnetic fieldstrengths using polarization by 6.7 GHz methanol masersaround the massive protostar Cepheus A HW2 and de-rive a line of sight (l.o.s) magnetic field strength ∼
23 mGat a distance of 300-500 AU from the central star. Themagnetic field strength is parametrized by the plasmabeta, β , which is the ratio of the thermal gas pressureto the magnetic pressure at the inner launching point l . Our simulations are so far limited to 1 < β < ∼
100 times weaker than esti-mated by conserving magnetic flux using observed valuesat 300 AU.Note that it is not only the field strength but also thefield distribution which affects the jet formation process,as it was shown by (Fendt 2006) by MHD jet forma-tion considering a wide parameter set of magnetic fieldand mass flux distribution along the jet launching area.As a result, simulations applying a concentrated mag-netic flux profile tend to be less collimated. We applyas initial field distribution the standard potential fieldsuggested by Ouyed & Pudritz (1997) ( § The stellar and disk luminosity
The massive young star produces a substantial lumi-nosity, which is supposed to dynamically change the out-flow structure. The dependence of the radiation force onthe stellar luminosity is parametrized by the Eddingtonratio Γ e , defined as the ratio of the radiation force dueto electron scattering to the central gravity (see Table1). The characteristic values of this dimensionless pa-rameter are obtained from the stellar evolution modelof Hosokawa & Omukai (2009). In order to study theimpact of radiation forces on the dynamics of outflows,we first launch a collimated jet from a disk wind viaMHD forces. Then, after the pure MHD jet has achieveda steady state, we initiate the radiation forces, and thencompare their impact on the MHD jet. Also, simulationsin which the radiation forces are considered from the be-ginning (which are computationally much more expen-sive) end up with a final structure of the outflow similarto the one obtained from the step-by-step method.We prescribe the radiation force by the line-drivingmechanism introduced by Castor, Abbott, & Klein(CAK). Such a force is parameterized by two physicalparameters k and α . The value of k is proportional tothe total number of lines. The quantity α can be con-sidered as a measure for the ratio of acceleration fromoptically thick lines to the total acceleration (Puls et al.2000). Depending on the selection of lines, for massiveOB stars typical values obtained for k range from 0.4-0.6, while α ranges between 0.3 and 0.7 (Abbott 1982).The process of line driving has also been applied to cata-clysmic variables (CVs) (Feldmeier & Shlosman 1999), aswell as to hot and luminous disks around Active Galac-tic Nuclei (AGN) (Proga et al. 2000, Proga & Kallman2004).Proga (2003) carried out numerical simulations driv-ing winds from hot luminous magnetized accretion disksof AGNs, assuming spherical geometry, an isothermalequation of state, and an initially vertical magnetic fieldstructure. Here we consider a potential field which ishour-glass shaped, and an adiabatic equation of state. Limitations of our model setup
Our paper, for the first time, provides a quantitativestudy of the interplay between radiative and MHD forceson outflows launched from the vicinity of young massivestars, applying high-resolution axisymmetric numericalsimulations. However, a few critical points can be raisedwhich may limit the applicability of our model and whichshould probably be considered in forthcoming investiga-tions.From the general point of view there is the lack of trueknowledge concerning a number of important parametersas discussed above. One important question is the loca-tion of the inner jet launching radius. If it is identicalwith an inner disk radius - where is that inner disk ra-dius located, if it exists at all? Is there a strong stellarmagnetic field which could open up a gap between thestellar surface and the disk as it is known for low-massyoung stars?Further, our disk model is taken as a boundary con-dition steady in time. As the star evolves, also the diskstructure and accretion rate may evolve in time. Thisquestion could only be answered by simulations solvingalso for the disk structure. We do not explicitly incor-porate heating and cooling but simply consider the adi-abatic expansion and compression in the jet.Another question is the existence of a stellar wind.We know that OB stars have strong mass loss in formof stellar winds during later stages of their life times.The derived mass loss rates are typically of the orderof 10 − M (cid:12) yr − . These winds are primarily radiationdriven via the line-driving mechanism (e.g. Kudritzki& Puls 2000, Owocki 2009). The velocities derived arehigh, ranging from ∼ − − and are usuallysupersonic. However, in case of high-mass young stars,no indication for such winds has been found so far possi-bly due to high obscuration. Nevertheless, a future studyshould implement the physical effect of a central stellarwind. BASIC EQUATIONS
For our study, we carry out axisymmetric numericalideal MHD simulations using the PLUTO code (Mignoneet al. 2007). We have modified the original code to in-corporate source terms treating the line-driven forcesfrom central star and disk, taking into account self-consistently the density and velocity distribution of theoutflow.The MHD code considers the following set of equations.That is the conservation of the mass, momentum, andenergy, ∂ρ∂t + ( (cid:126)v · ∇ ) ρ + ρ ∇ · (cid:126)v = 0 , (1) ρ ( ∂(cid:126)v∂t + ( (cid:126)v · ∇ ) (cid:126)v ) = −∇ P + 14 π ( ∇ × (cid:126)B ) × (cid:126)B − ρ ∇ Φ + ρ (cid:126)F rad , (2) ∂∂t ( ρE ) + ∇ · (cid:20) ρE(cid:126)v + ( P + B π ) (cid:126)v (cid:21) − (cid:126)B ( (cid:126)v · (cid:126)B ) = ρ (cid:104) −∇ Φ + (cid:126)F rad (cid:105) · (cid:126)v , (3)where ρ is the gas density, (cid:126)v the velocity vector, P thegas pressure, and (cid:126)B the magnetic field vector with thepoloidal and toroidal components (cid:126)B p , B φ . In order to in-clude the radiative forces (cid:126)F rad , the relevant source termshave been added in the momentum and energy conser-vation equation. The total energy density of the flow E comprises contributions from the internal energy (cid:15) , themechanical energy, and the magnetic energy, E = (cid:15) + v B πρ . (4)The gas pressure in the flow is related to the density as-suming an adiabatic equation of state with the adiabaticindex γ , P = ( γ − ρ(cid:15). (5)The evolution of the magnetic field is governed by theinduction equation, ∂ (cid:126)B∂t = ∇ × (cid:16) (cid:126)v × (cid:126)B (cid:17) . (6)We treat the ideal MHD equations without consideringresistive terms.In addition to the above set of equations the code obeysthe condition of divergence-free magnetic fields, ∇ · (cid:126)B =0, using the constraint transport method. Prescription of radiation forces
We do not explicitly consider radiative transfer, how-ever, we study the effects of momentum transfer by ra-diative forces on the outflow matter which is launchedby MHD processes from the underlying disk. The totalradiative force (cid:126)F rad comprises of four contributions - theacceleration due to continuum radiation from star (cid:126)f cont , ∗ and disk (cid:126)f cont , disk , respectively, and, similarly, due tospectral lines from star (cid:126)f line , ∗ and disk (cid:126)f line , disk , respec-tively, (cid:126)F rad = (cid:126)f cont , ∗ + (cid:126)f cont , disk + (cid:126)f line , ∗ + (cid:126)f line , disk . (7)In case of young massive stars, the stellar radiation issub-Eddington. This immediately implies that the con-tinuum force do not contribute in modifying the dynam-ics of the MHD outflow as its strength is much belowthe gravitational pull of the central star. However, lineforces could prove to be efficient in substantially enhanc-ing the continuum force by a so-called force multiplier,which is subject to complex theoretical studies of radia-tive transfer. This theory was developed first by Castoret al. (1975) who solved the radiative transfer equationsHD jets and radiation pressure 5from spectral lines in moving atmospheres. Accordingto their studies, the force due to line driving can be ex-pressed as a product of force due to continuum radiationand a force multiplier.The force multiplier depends on two parameters - k and α . A ”general” parametrization for the force multiplierwhich is independent of ion thermal velocities v th (seeEq A2 in the Appendix) has been introduced by Gayley(1995). In his formulation, the parameter k , initiallyintroduced by Castor et al. (1975), has been replaced bya parameter ¯ Q and the force multiplier can be expressedas follows, M ( T ) = k T − α = (cid:20) ¯ Q − α − α (cid:18) | ˆ n · ∇ (ˆ n · (cid:126)v ) | σ e cρ (cid:19) α (cid:21) , (8)where T is the optical depth parameter that depends onthe l.o.s velocity gradients (see eq. (A2)). The ”strongestform” of this parameterization requires the ansatz ¯ Q = Q , where Q being the line strength of the strongestline.Typically for evolved massive stars, the parameter ¯ Q lies in the range of 1000-2000 (Gayley 1995). For ourpresent study, we have applied the force multiplier con-sidering this ”general” parameterization in its ”strongestform”. Thus, we have the two parameters Q and α which define the force multiplier (Eq. 8) and hence theline forces. For simplicity, we assume constant values forthese line force parameters for a particular simulationrun. One of the important properties of the force mul-tiplier M ( T ) is that its values saturates to a maximumvalue M max when the medium becomes optically thinand the optical depth parameter approaches a minimumvalue, T ≤ T min . The typical value for the maximumforce multiplier is relatively constant, M max ∼ , de-pending only on the metallicity (Gayley 1995).The line force due to the central star is a product ofthe force due to continuum from the star and the forcemultiplier, (cid:126)f line , ∗ = (cid:126)f cont , ∗ M ( T ) . (9)The most difficult part for the numerical realization ofthe line forces is to calculate the proper l.o.s velocitygradients | ˆ n · ∇ (ˆ n · (cid:126)v ) | that appears in the formulationof force multiplier (Eq 8). Following the definition ofthis term by Rybicki & Hummer (1978), we express thisgradient in terms of the rate-of-strain tensor e ij ,ˆ n ·∇ (ˆ n · (cid:126)v ) = (cid:88) i,j e ij n i n j = (cid:88) i,j (cid:18) ∂v i ∂r j + ∂v j ∂r i (cid:19) n i n j , (10)where in general v i , r i , and n i are the components of ve-locity (cid:126)v , distance (cid:126)r , and the unit vector ˆ n , respectively.The different components of this tensor in cylindrical co-ordinates were calculated from Batchelor (1967). Progaet al. (1998) approximated the above sum to be equal tothe most dominant term, corresponding to the radial gra-dient of flow velocity along the spherical radius. Simula-tions with a more accurate algorithm of including all theother terms to determine the l.o.s. velocity gradient us-ing the above equation have also been carried out (Progaet al. 1999). These simulations show that the qualitativefeatures of winds are not changed as compared to the ap-proximate calculation of the gradients. Further, they are numerically very expensive. We can approximate thatthe region of consideration is far away from the star asthe central object is a point source. Therefore, the gra-dient can simply set to be equivalent to dV R /dR , where R is the spherical radius and V R is the gas velocity alongthe radius R . However, when using cylindrical coordi-nates instead of spherical coordinates we have to re-writeequation (10) by transforming R to r , and thereby addingfurther terms,ˆ n · ∇ (ˆ n · (cid:126)v ) = 11 + x (cid:18) ∂v r ∂r + x (cid:18) ∂v r ∂z + ∂v z ∂r (cid:19) + x ∂v z ∂z (cid:19) , (11)where x = z/r .We have applied two of the line force components indi-vidually as source terms in order to disentangle their ef-fects on the dynamics of a (pure) MHD jet launched fromthe underlying disk. These two components are from thecentral star alone and other due to underlying hot ac-cretion disk. Disks around massive young stars accretevery rapidly with rates of 10 − M (cid:12) yr − to 10 − M (cid:12) yr − .Such high accretion rates imply very high accretion lumi-nosities, in particular in their very inner regions. In or-der to calculate the line forces due to the disk luminosity,proper geometric factors have to be taken into account.We apply a formulation similar to Pereyra et al. (2000)whose details are given in the appendix. NUMERICAL SETUP
Grid setup & physical scaling
We perform axisymmetric ideal MHD simulations ofjet formation in the presence of radiative forces. Oursimulations are carried out on a grid of physical size ( r × z ) = (52 l × l ), where l denotes the physical lengthscale. The grid is divided into (512 × r < l and z < l the grid is uniform with a resolution of δr = δz =0 .
05, while for l < r < l or l < z < l thegrid is stretched with a ratio of 1.002739 and 1.001915 inradial or vertical direction, respectively. The remaining,very outer part of the grid is again uniform with cell size δr = 0 .
125 and δz = 0 .
20, respectively.Pure MHD simulations would be scale-free. However,when radiation is considered a physical scaling of the dy-namical variables becomes essential. Three scaling pa-rameters in physical units are used – the length scale l ,the base flow density ρ and the Keplerian velocity v ,at the inner launching point (see § r c = r cgs l , z c = z cgs l , ρ c = ρ cgs ρ , v c = v cgs v ,p c = p cgs ρ v , B c = B cgs B = B cgs (cid:112) πρ v . (12)The force multiplier defined in equation (8) can be ex-pressed in code units as follows M c ( T c ) = (cid:20) Q − α − α (cid:18) v σ e cρ l ρ c (cid:12)(cid:12)(cid:12)(cid:12) dv l dl (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) α (cid:21) . (13)Quantities with subscript ’c’ are obtained from simu-lations in the dimensionless form where as the quantitieswith subscript ’cgs’ are the ones required in the physi-cal units. We measure the time t c in units of l /v . Wedenote the number of rotations at the inner launchingpoint l as N rot = l / πv = t c / π Using these normalizations the conservation of momen-tum (Eq. 2) can be written in dimensionless code units, ρ c (cid:18) ∂(cid:126)v c ∂t c + ( (cid:126)v c · ∇ c ) (cid:126)v c (cid:19) = 2 β (( ∇ c × (cid:126)B c ) × (cid:126)B c ) −∇ c (cid:126)P c − ρ c ∇ c Φ c + ρ c ( (cid:126)F rad c ) , (14)where the plasma beta, β = (8 πρ v ) /B , specifies theinitial strength of the (poloidal) magnetic field at l andthe radiation force (cid:126)F radc = (cid:126)F radcgs v /l . (15)Figure 2 displays the numerical setup used for oursimulations. We also show the dynamically importantforces that a fluid parcel experiences in an outflow froma young massive star - the gravitational force (cid:126)F gravity by the central star, the Lorentz force components paral-lel (cid:126)F Lorentz , (cid:107) (accelerating) and perpendicular (cid:126)F Lorentz , ⊥ (collimating), the thermal pressure gradient (cid:126)F pressure ,and the centrifugal force (cid:126)F centrifugal , which plays a vi-tal role in particular for accelerating the flow from thedisk surface. The essential point of this paper is thatwe also consider radiation forces from the star and thedisk ( (cid:126)F radiation , eq. (7)). We consider the most domi-nant radiative source terms which are those due to theline driving mechanism from the central massive star andthe underlying (hot) disk (see Fig 1). Initial conditions
We model the launching of the wind from an accretiondisk representing the base of the outflow. The gravita-tional potential Φ is that of point star located at a slightoffset ( − r g , − z g ) from the origin to avoid singularity atr=z=0, Φ( r, z ) ∝ (cid:0) ( r + r g ) + ( z + z g ) (cid:1) − . , (16)where r g = z g = 0 .
21 in all our simulations.Initially we prescribe a hydrostatic equilibrium with adensity distribution ρ ( r, z ) ∝ (cid:0) ( r + r g ) + ( z + z g ) (cid:1) − . (17)(Ouyed & Pudritz 1997). The thermal pressure followsan polytropic equation of state P = Kρ γ , with γ = 5 / K = ( γ − /γ is determined by theinitial hydrostatic balance between gravity and pressuregradient. To begin with, all velocity components are van-ishing. Such an initial setup is typical for jet launchingsimulations (e.g. Ouyed & Pudritz 1997, Krasnopolskyet al. 1999, Fendt 2006). The initial hot hydrostatic dis-tribution of density and temperature does not play a sig-nificant role in determining the flow structure as it iseventually washed out by the cold [i.e. thermal pressureunimportant] MHD flow that is launched from the base.The initial magnetic field is purely poloidal (in the( r, z )-plane) with a distribution derived from the poten- Figure 2.
Schematic view of the numerical setup along with defi-nition of the boundaries conditions applied for most of the simula-tion runs. The lines indicate poloidal magnetic field lines anchoredin the underlying accretion disk which is in slightly sub-Keplerianrotation. The dark black spot represents a fluid element frozenonto that field line. Vectors indicate and denote the various forcesacting on the fluid element. See text for a detailed description ofthe different dynamically important forces. tial field (cid:126)B p = ∇ × A φ with A φ = (cid:112) r + ( z + z d ) − ( z d + z ) r (18)(Ouyed & Pudritz 1997), where z d is considered as thedimensionless disk thickness. Such a magnetic field con-figuration is force free. Together with the hydrostaticdensity distribution this implies an initial setup in forceequilibrium. The initial magnetic field strength is pre-scribed by choice of the plasma-beta β at the innerlaunching point. Boundary conditions
To choose the correct physical boundary conditions isof utmost importance for numerical simulations as theydescribe the astrophysical system under consideration.In the present setup, we have to deal with four boundaryregions (see Fig. 2).
Axial boundary
Along the jet axis an axisymmetric boundary conditionis applied. The normal and the toroidal components ofvector fields (i.e. B n , B φ , v n , v φ ) change sign across theboundary, whereas the axial components are continuous.HD jets and radiation pressure 7The density and the pressure are copied into the ghostzones from the domain. Equatorial Plane boundary
The ”disk boundary” is divided in two regions - theinner gap region extending from the axis to the innerlaunching radius, r < l , and further out the disk region, r ≥ l , from where the outflow is launched.The setup of this boundary is the most crucial for theoutflow simulation. This is an ”inflow” boundary withthe boundary values determining the inflow of gas andmagnetic flux from the disk surface into the outflow. Spe-cial care has to be taken to consider the causal interactionbetween the gas flow in the domain and in the ghost cells.In order to have a causally consistent boundary condi-tion (Bogovalov 1997, Krasnopolsky et al. 1999, Porth& Fendt 2010), we impose the four physical quanti-ties ρ, Ω F , P, E φ . The toroidal electric field vanishes, E φ = ( (cid:126)v × (cid:126)B ) φ = 0, as result of the ideal MHD ap-proximation. Thus, poloidal velocity and the poloidalmagnetic field are parallel at the boundary, (cid:126)v p || (cid:126)B p . Theangular velocity of the field line (Ferraro’s iso-rotationparameter) Ω F , which is one of the conserved quanti-ties along the field line is fixed in time along the bound-ary. We have chosen a Keplerian profile along the diskΩ F ( r ) ∝ r − . . The poloidal velocity at the boundaryis ”floating”, i.e. copied from the first grid cell into theghost cell each time step. Thus, the mass flux at the diskboundary is not fixed but consistently derived by causalinteraction between the outflowing gas and the boundaryvalue.The inner gap region is prescribed as hydrostatic pres-sure distribution. However, gas pressure gradient, grav-ity, and centrifugal acceleration is considered for the ra-dial force-balance in the disk, leading to a sub-Keplerianrotation v φ ( l , v kep = √ χ with χ <
1. Solving for the radial balance delivers thedensity profile along the boundary ρ disk ( r, z ) = (cid:18) − χ (cid:19) ρ ( r, z ) . Therefore, there is a density contrast between disk andinitial corona (the domain). In order to avoid numericalproblems at the interface of between gap and ”inner disk”(ie. the inner jet launching radius), we smoothen thistransition using the Fermi function, considering χ = (cid:18) χ − r − (cid:19) . The Fermi function is resolved by 16 grid cells. Thepressure distribution along the boundary is fixed so asto have a cool disk which is three times denser than theinitial corona. Thus, the wind emerging from such acool disk has a similarly cool temperature (lower thanthe initial corona). Quantitatively, the average initialcoronal temperature within the inner jet launching areais around 5 × K, while the temperature of the disksurface (i.e., the boundary) in the same area reduces to5 × K at 10 AU.
Table 1
Summary of the dimensionless parameters to study the impact ofradiation forces on the outflow dynamicsParameter DefinitionEddington ratio Γ e ≡ σ e L ∗ πcGM ∗ (proxy for stellar luminosity)Luminosity ratio µ ≡ L acc L ∗ = GM ∗ ˙ M acc R ∗ L ∗ (proxy for disk luminosity)Radius ratio Λ ≡ R ∗ l (proxy for stellar radius)Initial plasma beta β ≡ πP B (proxy for magnetic field strength)Line-force parameters Q and α prescribes M ( T ) using Eq. (8) The magnetic field in the boundary is allowed to evolvein time to avoid any current sheets along the interface.When the simulation starts, the rotating disk winds-upthe poloidal field and induces a toroidal field component.Constraining the field line angular velocity Ω F to be con-stant in time, we need to adjust the rotational velocityof the gas in the boundary, v φ = r Ω F + ηρ B φ , (19)where the mass load of a field line η = ( ρv p ) /B p is againa conserved quantity in stationary MHD. Outflow boundaries
The right and top boundaries (see Fig. 2) are definedas outflow boundaries. The canonical outflow conditions(zero gradient across the boundary) are imposed to allscalar quantities and vector components, except for themagnetic field. The toroidal magnetic field component B φ in the boundary is obtained by requiring the gra-dient of the total electric current I to vanish, whereasthe poloidal field components are estimated by settingthe toroidal current density j φ to zero (Porth & Fendt2010). Having this new outflow condition implementedin the standard PLUTO code, we ensure that there is noartificial collimation due to boundary effects. PARAMETER SURVEY
Choice of governing parameters
The three main parameters we apply for a compre-hensive study are the (i) central stellar mass M ∗ , (ii)plasma-beta β , and (iii) the line force parameter α (seeTable 1). We also find that the magnitude of the den-sity at the launching base of the flow ρ , which is a freeparameter in our simulations, strongly affects the flowcharacteristics (see discussion in § (cid:12) to 60 M (cid:12) . The size and Figure 3.
Variation of the dimensionless parameter Γ e with stel-lar mass M ∗ from Hosokawa & Omukai (2009). The solid black line indicates the stellar mass evolution in the Hosokawa model. Stars represent the Eddington ratio Γ e for different stellar masses in ourparameter survey § vertical dashed lines mark differentevolutionary stages for massive young star accreting rapidly at therate of 10 − M (cid:12) yr − . The curve is obtained using tabulated evolu-tionary parameters kindly provided by Takashi Hosokawa (privatecommunication). the luminosity of these stars are obtained from the lit-erature stellar evolution models (Hosokawa & Omukai2009). The dimensionless parameter that controls theluminosity of the central star is Γ e (see Table 1). Itsvariation with the stellar masses considered in our modelis shown in Fig. 3. (Also see § β = 1 . , . , .
0. We fixthe density at the base of the flow to the same value asfor the simulations for different stellar mass. Thus, for30 M (cid:12) star and an inner jet launching radius of 1 AU, thedifferent β values correspond to a poloidal magnetic fieldof 11.5 G, 6.6 G, or 5.1 G at the inner launching point, re-spectively. The detailed description of the different sim-ulations is shown in Table 2. The jet is launched bypure MHD processes. As it achieves a steady state af-ter say N MHD inner rotations, we switch on the radiativeline forces by the central star and run the simulation foranother N RAD inner rotations.The stellar line force in our model is quantified by theparameters Q and α . In principle, these parameters de-pend on the degree of ionization in the flow, however, forsimplicity, we neglect this dependence and assume thattheir values remain the same for a particular simulationrun. The line force is independent of Q , however, theits dependence on α is strong Gayley (1995). We fix theline strength parameter Q = 1400 for all runs. We carryout simulations with three values of α = 0 . , . , . α seems to be small, theylead to considerable variation in the magnitude of theline forces (Gayley 1995). These values of α and Q areconsistent with the empirical values obtained for evolved massive stars (Abbott 1982,Gayley 1995). Since, we haveno empirical values for these parameters during the for-mation stage, the present model uses similar values ofline force parameters as obtained from evolved massivestars. Quantifying the degree of collimation
There are in principal several options to quantify thecollimation degree of jets. Fendt & ˇCemelji´c (2002) sug-gested to compare the mass fluxes in radial and verti-cal direction of the jet as a measure of the jet collima-tion. However, in general the choice of a ”floating” inflowboundary condition and the subsequent time-dependentchange of mass fluxes makes it difficult to use them asa measure for the degree of collimation. Therefore, inthis work, we quantify the jet collimation by the openingangle φ of the field lines (which are equivalent to velocitystreamlines in a steady state). For comparison, we mea-sure the opening angle along a certain field line at twocritical points along that field line viz. the Alfv´ e n point φ A and the fast magneto-sonic point φ F . Note that thisis not the asymptotic jet opening angle.In order to estimate the amount of de-collimation∆ ζ [%] by radiation forces, we measure the angular sep-aration of field lines resulting from N MHD disk rotationsfrom those with after N RAD disk rotations,∆ ζ [%]( s ) = 100 (cid:18) φ ( s, N RAD ) − φ ( s, N MHD ) φ ( s, N MHD ) (cid:19) , (20)where s measures of path along the field line. A positivevalue of ∆ ζ [%] implies de-collimation whereas negativevalues would imply collimation of the jet by radiativeforces. The outflow density at the disk surface - ρ The radiative force term F radc is a product of the forcemultiplier and the continuum force. The force multiplier(Eq. 8) depends on the physical scalings v , l and ρ (see Eq. 13). This does imply that each simulation runis unique to the chosen set of scaling parameters for aparticular type of star.In case of massive stars these values are currently dif-ficult to measure very close to the star (see § − to 10 − M (cid:12) yr − .The density ρ at the base of the flow ( z ∼
0) can beestimated by the mass flux launched from the base perunit time, ˙M out = 2 πρ v l (cid:90) r c , max r c , min r / − q c dr c , (21)where the density at the base of the flow is ρ (r c , z c = 0) = ρ r − qc . From Eq. (17), q = 3/2. We assume the matteris launched from the disk surface with Keplerian speed(i.e. v z ( r c , z c = 0) = v r − / ). The physical scalingfor the density, velocity and lengths are ρ , v and l HD jets and radiation pressure 9
Table 2
Parameter study of stellar radiation line-force effects on MHD disk jets for different field strength and different stellar mass. We apply thesame physical density at the jet base ρ = 5 . × − g cm − , the same inner launching radius l = 1 AU, and the same line-forceparameter Q = 1400 for all runs. The simulations are performed for N MHD = 319 inner rotations in pure MHD, followed by and N RAD = 319 rotations, with switched-on radiative forces. A measure for the jet opening angle is given at the MHD critical points, φ A , φ F and along the field line rooted at r = 5 . AU . ∆ ζ [%] denotes the percentage difference of the opening angles between the steady-stateflows for pure MHD and including radiation force along the field line rooted at the same radius. The vertical mass flux ˙ M vert [M (cid:12) yr − ] ismeasured along the top cells in the domain. The first ten simulation runs in the table apply a ”floating” boundary condition for theinjection velocity, while for last three runs a fixed-mass-flux boundary condition is applied (see § M ∗ [M (cid:12) ] β α φ A φ F max[∆ ζ [%]] ˙ M vert [M (cid:12) yr − ]M30a055b1 30 1.0 0.55 21.98 16.45 5.0 3 . × − M30a055b3 30 3.0 0.55 25.63 20.68 17.0 2 . × − M30a055b5 30 5.0 0.55 26.05 23.40 34.5 2 . × − M30a055b5-c 30 5.0 0.55 25.55 23.01 33.9 2 . × − M20a055b5 20 5.0 0.55 20.44 15.57 5.8 1 . × − M25a055b5 25 5.0 0.55 23.27 19.94 21.9 2 . × − M45a055b5 45 5.0 0.55 26.84 24.35 37.7 3 . × − M50a055b5 50 5.0 0.55 28.86 26.03 42.0 3 . × − M60a055b5 60 5.0 0.55 31.45 29.25 56.0 3 . × − M60a060b5 60 5.0 0.60 23.73 21.48 28.7 3 . × − M60a065b5 60 5.0 0.65 20.93 17.28 11.7 2 . × − M25a055b3inj 25 3.0 0.55 45.84 33.87 30.2 1 . × − M30a055b3inj 30 3.0 0.55 47.71 36.97 37.2 1 . × − M50a055b3inj 50 3.0 0.55 48.76 39.22 42.1 1 . × − respectively, while r c is the non-dimensional length unit.Thus, the density scaling ρ can be expressed as ρ = ˙M out πv l (cid:20) ln (cid:18) r c , max r c , min (cid:19)(cid:21) − . (22)Using typical observed mass outflow rates for young mas-sive stellar jets, we calculate ρ ∼ − − − g cm − .The physical value of the density at the base of the flowis in the denominator of Eq. 13, and therefore affects themagnitude of the line-driving force significantly. Increas-ing the density by, say, a factor of 10, decreases the lineforce at least by a factor of 10 α . For the values of α considered her, this magnification could be as large as 3.Such a change in the radiative force has clearly a notableimpact on the outflow dynamics, in particular on its de-gree of collimation. The description of our simulationswith different base density is given in Table 3. RESULTS AND DISCUSSION
The radiation field in massive star forming regions mayplay a crucial role in modifying the dynamics of outflowsand jets. In order to disentangle the effects of radiativeforces from the pure MHD jet formation, we decided tofollow a two-step approach. We first i) launch a pureMHD disk jet and wait until it has reach a steady state(after about N MHD = 320 rotations). We then ii) switchon the radiation line-forces and allow the jet flow to fur-ther develop into a new dynamic state. In some cases anew steady-state can be reached, in other cases an un-steady solution develops.For the radiation forces we consider line-forces exertedby i) the luminous massive young star, and those by ii)the surrounding accretion disk.
Jet de-collimation by the stellar radiation field
The radiation line-force from the central point star isdetermined by applying Eq. (9). For a detailed physicalanalysis of the effects of line-forces from the star alone,we have chosen simulation M30a055b5 as reference sim-ulation run (Table 2). The reference run is parametrizedfor a stellar mass of 30M (cid:12) , surrounded by an accretiondisk with an inner jet launching radius of 1 AU, and adensity of 5 . × − g cm − at this radius. The initialpoloidal magnetic field strength is fixed by a plasma- β = 5 .
0. The radiative forces from the star are definedby the line force parameters Q = 1400 . α = 0 . M ( T )as well as the specific stellar line force for the referencesimulation are both shown in Fig. 4. The magnitude ofthe force multiplier peaks in the top-right low-density re-gions of the flow. This is expected as the force multiplierincreases with decreasing density as shown in Eq.9. Inorder to calculate the true radiation force, the force mul-tiplier must be convolved with the continuum radiationfrom the star, resulting in a line-force distribution witha maximum close to the central star (see Fig. 4).The time evolution of the emerging jet † is shown inFig. 5, where we display the vertical jet velocity v z ( r, z )and the poloidal magnetic field structure for referencesimulation M30a055b5. We clearly see the change fromthe pure MHD flow ( top panels) to the situation whenstellar radiation forces are considered ( bottom panels).Material which is injected from the disk surface ( bottomboundary ) is ”frozen” on these field lines (ideal MHD).Initially, the plasma is accelerated magneto-centrifugallyand gains substantial speed, producing a bow shock asit propagates. The bow shock leaves the domain after ∼
60 rotations. Inertial forces of the outflowing massflux induce a strong toroidal magnetic field component † Full movies are available for download at Table 3
Parameter study of stellar radiation line-force effects on MHD disk jets for different jet density at the base of the flow with a stellar massof 30 M (cid:12) . For all runs, plasma- β = 5, l = 1 AU and line-force parameters Q = 1400 and α = 0.55.Runs ρ [g cm − ] φ A φ F max[∆ ζ [%]] ˙ M vert [M (cid:12) yr − ]M30a055b5d1 3 . × − . × − M30a055b5 5 . × − . × − M30a055b5d2 1 . × − . × − M30a055b5d3 5 . × − . × − Table 4
Parameter study of disk radiation line-force effects on MHD disk jets. All simulation runs apply the same stellar parameters andline-force parameters as M30a055b5, however a different inner launching radius of only 0.1 AU. Additionally, two more dimensionlessparameters for the disk radiation force µ = 0 . . ρ [g cm − ] max[∆ ζ [%]] ˙ M vert [M (cid:12) yr − ] CommentsDisk1 5 . × − -1.8 6 . × − reaches a steady stateDisk2 5 . × − . × − remains unsteady Figure 4.
Contours of force multiplier M( T ) and the specific lineradiation force of a central star [in physical units] for the referencesimulation M30a055b5 after N rot = 638 inner rotations. resulting in magnetic hoop stresses which self-collimatethe magnetic field structure together with the hydrody-namic mass flux. Eventually, when all the dynamicalforces along the field line are balanced again, the flowachieves a steady-state. The steady state magnetic fieldconfiguration obtained after the pure MHD flow is shownas white dashed lines in Fig. 5.At this stage we ”switch on” the line-forces of the cen-tral star. Immediately, the emergence of a fast axial flowis visible. This flow is unsteady, first forming a knottystructure, and then, when approaching the upper bound-ary, stabilizes to a steady axial flow. Jet de-collimationby stellar radiation forces is indicated by the fact thatthe red solid field lines in the bottom panels of Fig. 5 doopen-up significantly as compared to field lines in steady- state pure MHD simulation.When the radiative force is switched on, a shock frontbegins to propagate into the steady MHD jet (Fig. 5).This happens, because the additional radiation forceslead to an initial disturbance at the base of the flow,which is then propagated outwards. The effect of line-forces resulting in a series of propagating shocks is bestseen in the poloidal velocity profile along a magnetic fieldline (see Fig. 6). The series of shocks eventually prop-agate out of the domain, and the flow attains a steadystate again. However, the asymptotic outflow velocitywhich is then achieved is enhanced by a factor 1.5 - 2compared to the pure MHD flow.We have also carried out a run, M30a055b5-c, that in-cludes the radiation force right from the beginning alongwith the MHD flow. This flow evolves into the sameconfiguration as our reference run where radiation forceis ”switched on” later. This proves that the initial con-ditions do not affect the final state of the flow. In ourstudy, we prefer to use the step-by-step approach as it iscomputationally less expensive.In summary, we find that radiation forces modify theMHD disk jet essentially in terms of collimation, but alsoin terms of acceleration and terminal speed. Analysis of the force-balance in the outflow
Here we investigate the main forces affecting the out-flow dynamics, comparing the magnitude of various forceterms calculated at each grid point along the field linerooted at r = 5.0 AU for simulation run M60a055b5. Fig-ure 7 shows such a comparison of specific forces pro-jected parallel and perpendicular to the field line beforeand after considering line-forces due to stellar radiation.With respect to acceleration, the most striking featureis the enhanced specific pressure gradient force. Whenradiative forces are considered, the steady MHD flow ma-terial at higher altitudes is adiabatically compressed fromthe underlying accelerated material leading into a shockthat propagates out of the domain. The resultant flowhas higher temperature or an increase in thermal pres-sure. We also find that the gradient of the thermal pres-sure higher above in the flow may increase substantially,HD jets and radiation pressure 11
Figure 5.
Evolutionary sequence of the reference simulation M30a055b5. Shown is the jet vertical velocity distribution (color coding)and the poloidal magnetic field lines. The color bar represents the velocity scale in km s − for each row of panels. The series of imagesare taken at subsequent inner rotations as mentioned on the top of each image. The solid red lines are the poloidal magnetic field lines.To illustrate the impact of radiation forces, we also show the field lines obtained after the steady-state MHD flow for comparison as whitedashed lines in the lower four panels. such that thermal pressure force and Lorentz force maybecome comparable. In terms of collimation, it is evidentfrom Fig. 7 that the profile of the perpendicular compo-nent of the pressure gradient force along the field line be-comes flatter than in the pure MHD flow case. Thus, thepressure force which was not important for pure MHDflows, now becomes a significant factor in governing thedynamics of the initial flow acceleration (disk wind).Contrary, when the line-driving force is switched-onin the MHD jet, the centrifugal force becomes reduced by an order of magnitude particularly at high altitudesin the flow. In order to comprehend this variation ofcentrifugal force, it is useful to compare the conservedquantities of MHD, in particular the angular velocity ofthe field line Ω F . In the pure MHD simulation Ω F isconserved throughout the simulation. This remains truealso with radiative forces included, simply because Ω F is fixed as boundary value (see § Figure 6.
Poloidal velocity along the field line rooted at r = 3AU for the two simulations M30a055b5 top panel and Disk2 bottompanel . The distance along the field line s is measured in AU. Thevarious colored solid curves indicate the velocity profile at differentinner disk rotations N rot ( blue for N rot = 325, green for N rot = 333, red for N rot =340, cyan for N rot = 510, and magenta for N rot =638). The black dotted line in both panels correspond to the lasttime step for the pure MHD flow. The solid yellow line shown inthe bottom panel indicates the poloidal velocity profile after 1500inner rotations. accelerate the flow (close to the base) to higher poloidalvelocities. The magnetic flux at the base is fixed as aboundary condition. Hence, the ratio of mass load todensity η/ρ in Eq. (19) increases, and since the toroidalmagnetic field B φ is negative, also the azimuthal veloc-ity decreases. This eventually leads to a decrease of thespecific centrifugal force when line-forces are considered.The decrease in centrifugal force further implies that theoutflow rotates slower, and that the acceleration close tothe base is not only controlled by magneto-centrifugalforces, but also by the thermal pressure force and theradiation force (see Fig. 7).Close to the base of the flow, Lorentz forces are dy-namically not important. However, they peak at theAlfv´ e n point of that particular field line. Beyond theAlfv´ e n point, the Lorentz force becomes important andcompetes with other forces to govern the dynamics of theflow. When radiative force are considered, the Lorentzforce close to the base of the flow and also at higher alti-tudes do not show significant deviations from its valuesin the pure MHD case.Line-driven radiative forces have considerable impacton the poloidal electric current distribution and alsoon the critical MHD surfaces in the jet. The contourplots shown in Fig. 8 compare the total poloidal current I = rB φ = (cid:82) (cid:126)j p d (cid:126)A and position of critical surfaces be-fore and after considering line-driven forces in the refer-ence simulation M30a055b5. The figure shows that byadding radiative forces (stellar luminosity), the corre-sponding poloidal electric current density contours areshifted closer to the base of the flow, implying a lowertoroidal field strength B φ in these jets. We understandthis results as a consequence of lack of jet rotation, thus less effective induction of toroidal magnetic field.Jet collimation in the conventional Blandford-Paynepicture is caused by magnetic hoop stresses (- B φ /r ) ofthe azimutal magnetic field. Whereas, the larger toroidalmagnetic pressure gradient aids in de-collimating theoutflow. A lower B φ would not only weaken the hoopstress but also reduce the magnetic pressure ( ∼ B φ ).The decrease of toroidal magnetic pressure would lead toa lower magnetic pressure gradient force. So, a balancearises between the hoop stress that collimates the flowand the toroidal pressure gradient that de-collimates it.We observe that the field lines de-collimate to a widerconfiguration on addition of radiative forces, implyingthat the decrease of toroidal magnetic field has more im-pact on reducing of hoop stresses than reducing the mag-netic pressure gradient, thus eventually de-collimatingthe flow.Considering radiative forces also results in lower-ing the location of MHD critical surfaces (viz. slowmagneto-sonic, fast magneto-sonic and Alfv´ e n surface)(see Fig. 8). The pure MHD flow is launched marginallysuper-slow and sub-Alfv´ e nic. The flow near the axis is,however, sub-slow, due to the boundary condition of con-serving the initial hydrostatic density and pressure dis-tribution in the gap region between axis and disk. Theflow speed at the critical points depends on the magneticflux and mass density in the flow. These quantities re-main approximately unaltered when radiative forces areconsidered, thus the magneto-sonic wave speed remainssimilar as well. Since now the radiative line-forces ac-celerate the wind further to higher velocities, the criticalsurfaces shift to lower altitudes in the flow (i.e. close tothe base of the flow).We finally note that radiative forces have both a directas well as an indirect effect in modifying the collimationproperties of the flow. The radiation force from the star directly affects the flow dynamics at its base and close tothe star simply by transfer of momentum (Eq. 2), even-tually leading to a flow de-collimation. Indirectly , theradiation field also enhances the thermal pressure in theflow (Eq. 3) such that the specific pressure force becomescomparable to the Lorentz force - further de-collimatingthe flow. The resulting opening angles as measured fromour reference run at various critical points along the fieldlines are listed in Table 2. For M30a055b5, the fi-nal state of the jet has a steady mass outflow rate of2 . × − M (cid:12) yr − . The opening angles are 26 ◦ and 23 ◦ at Alfv´ e n and fast point respectively for the field linewhich has a foot point r = 5 AU. The maximum per-centage separation ∆ ζ [%] between the opening angles ofthis field line in pure MHD flow and that for the flow withradiative forces is 34.5%. A significantly positive valueof ∆ ζ [%] indicates flow de-collimation. In the followingsections, we describe the effects of various physical pa-rameters that could affect the dynamics of the outflow. Radiation field and magnetic field strength
Here we discuss the interrelation between radiativeforces and a variation of magnetic flux and their com-bined effect on outflow collimation. We compare sim-ulation runs M30a055b1, M30a055b3, and M30a055b5,all assuming the same stellar mass 30M (cid:12) , an inner jetlaunching radius of 1 AU, and a density at the base ofHD jets and radiation pressure 13
Figure 7.
Comparison of all specific force terms projected parallel [left] and perpendicular [right] to the field line rooted at r=5.0 AU forthe run M60a055b5. Colors specify different specific force terms in dyne cm g − - the Lorentz force ( green ), the gas pressure gradient( red ), and the centrifugal force ( blue ). The dashed lines corresponds to the respective force terms for the pure MHD run M60a055b5 (atsteady-state), while the solid lines represents the forces for the final time step of the simulation including radiative forces. The brown solidline represents stellar radiation force term for the final time step. The black and magenta vertical solid lines mark the position of theAlfv´ e n and the fast magneto-sonic surface for the steady state flow including radiative forces, whereas the corresponding dashed lines arefor the pure MHD flow. the flow of 5 . × − g cm − , while the initial magneticfield strength parametrized by β ranges from 5.1 G to11.5 G.The (initial) pure MHD runs exhibit characteristic dif-ferences. A lower β provides faster jets which are morecollimated, simply because the larger field strength pro-vides larger Lorentz forces to collimate and also accel-erate the flow more efficiently. A lower field strengthimplies a lower magneto-sonic wave speed, resulting inmagneto-sonic surfaces located closer to the disk surface.Thus, we cannot simply compare the opening angles atthe critical surfaces as their positions vary in the pureMHD runs with different field strengths. Instead, wequantify the actual change of opening angle with andwithout considering radiative force ∆ ζ [%] (i.e. in per-centage), as profile along the field line.This measure is shown in Fig. 9 for a field line rootedat r = 5 AU. We see that ∆ ζ [%] is significantly posi-tive, indicating that radiative forces do considerably de-collimate the flow. However, the effect of de-collimationis enhanced when the magnetic field strength is reduced(i.e. an increase of β ). The asymptotic value of ∆ ζ [%]along a field line changes from 35% to 5% when the initialmagnetic field strength is increased by a factor of two.The above analysis for different magnetic field strengthsuggests that in the strong field case ( β = 1) jet colli- mation is controlled by the magnetic forces alone. How-ever, when we decrease the field strength (from 11.5 G to5.1 G as normalized at 1 AU) the stellar radiation line-force begins to compete with the magnetic forces. Theresulting jet has a much wider field structure and outflowopening angle as it is de-collimated as compared to itspure MHD counterpart (see Fig. 5). We conclude thatradiative forces from the central star (for the chosen pa-rameters) will dominate the magnetic effects in control-ling the flow collimation for field strengths (cid:46) § Impact of the line-force parameter α The stellar line-force is a non-linear function of α . Gay-ley (1995) showed that a significant change in the massflux can be obtained with very slight change in α , indicat-ing that this parameter is very sensitive for calculationof radiative line-forces.Simulation runs M60a055b5, M60a060b5, andM60a065b5 apply the same magnetic field strength( β = 5 . α = 0 . − .
65. A first result is that the increase of4
Figure 8.
Contours of the total poloidal electric current distribu-tion ( left ), and the location of the MHD critical surfaces ( right ).In both panels, the pure MHD flow is represented by black solidlines , while the radiative MHD flow is shown with red solid lines . Figure 9.
Flow collimation along the field line rooted at r=5AU. Shown is the profile of the percentage separation of field linesbetween the pure MHD flow and the radiative MHD flow [∆ ζ [%]]as measure of flow collimation. The solid , dashed and dot dashed lines are for three different plasma- β = 5 . , . , .
0, respectively. α from 0.55 to 0.60 decreases the mass outflow ratesignificantly by ∼ § α become more de-collimated thanthe others with higher α .The above trend indicates that the efficiency of the line Figure 10.
Jet collimation, jet density and radiation force.Shown is the opening angle (i.e. the field inclination) at the Alfv´ e npoint ( stars ) and the fast point ( dots ) of the field line rooted atr=5 AU for runs with different α ( top panel ), and different jet basedensity ρ ( bottom panel ). The dotted and solid lines correspondsto the opening angle at the Alfv´ e n point, and the fast point of thesame field line after the pure MHD flow with an initial plasma- β = 5.0. driving mechanism is increased with lower α . Physicallya lower α implies a higher contribution of optically thinlines in accelerating the flow. Thus, the dominance of lesssaturated (i.e. less self-shadowed) lines in acceleratingthe flow results in an efficient line driving.In summary, we find that the runs with lower α havehigher mass flux and are less collimated. Our resultsfrom studying different α are in qualitative agreementwith the analytical results from Gayley (1995). Impact of the density ρ In this section, we describe how collimation and ac-celeration in the outflow are altered with the density atthe base of the flow. As mentioned in § ρ = 3 × − g cm − For this run,the opening angles measured at critical MHD surfaces arelarger compared to our reference run. Quantitatively, theopening angles at the Alfv´ e n and the fast surface are 30 ◦ and 28 ◦ , respectively. Also, the maximum percentagechange in field line opening angle ∆ ζ [%] increases from35% for the reference run to a staggering high value of48% for field line rooted at r = 5 AU. In case of a highjet base density ρ , the change in field line opening angle∆ ζ [%] reduces to 4%, with the opening angles at theAlfv´ e n and fast surface being 20 ◦ and 15 ◦ , respectively(Fig. 10).These results correspond to an inverse correlation ofthe density at the base of the flow with the force multi-plier M ( T ) (see Eq. 13). Higher densities result in an op-HD jets and radiation pressure 15tically thick environment, thus increasing the the opticaldepth parameter T , and reducing the force multiplier.Thus, the radiative line-driving force approaches thelimit corresponding to the continuum radiation force. Es-sentially, the continuum force for typical massive youngstars is weak compared to other dynamically importantforces (for e.g. gravity).In summary, we observe the outflow density as one ofthe leading parameters to govern the dynamics of theoutflow by radiation forces. However, observationallythis density is very difficult, even impossible to deter-mine. We thus have to rely on certain assumptions,or may apply estimates of the outflow flux to calculatethis density. Our parameter survey considers density val-ues which are consistent with the observed mass fluxes.From our studies, we observe in general that for densities ρ < − g cm − the line-driving force from the centralstar is very efficient in accelerating and de-collimatingthe outflow. However, a denser environment will dilutethe influence of the line-driving mechanism. Stellar mass evolution and outflow collimation
Motivated by the hypothesis of Beuther & Shepherd(2005) we have studied the outflow dynamics and col-limation for a sequence of different stellar masses (viz.from 20 M (cid:12) to 60 M (cid:12) ). The change in stellar mass im-plies a change in the central luminosity. Fig. 11 showsthe variation of the opening angle at the MHD criticalpoints in simulation runs with central stellar mass. Itcan be seen that the opening angle varies from 20 ◦ to32 ◦ for stellar masses from 20 M (cid:12) to 60 M (cid:12) . The curverises linearly up 30 M (cid:12) , but then the opening angle doesnot change considerably for increasing mass.For the physical parameters of the stellar evolutionsuch as stellar luminosity and radius, we have appliedthe bloating star model of Hosokawa & Omukai (2009)assuming an accretion rate of 1 . × − M (cid:12) yr − (seefig. 3). In this model for massive young stars, the stel-lar luminosity increases rapidly from 6 M (cid:12) to 30 M (cid:12) . Astar with high accretion rate undergoes swelling and isconsidered to bloat up to 100 R (cid:12) . It is the entropy dis-tribution in these stars which causes them to expandin size that much. For stars <
10 M (cid:12) , a thin outerlayer absorbs most of the entropy from its deep interiors,thereby increasing the entropy in this layers. The rapidincrease of entropy in the outer layer of the star causesthe star to increase in size (Hosokawa & Omukai 2009).When the stellar mass reaches ∼
10 M (cid:12) , the star thenbegins to shrink (Kelvin-Helmholtz contraction) until itsmass reaches ∼
30 M (cid:12) . After that, the star approachesthe main sequence and then follows the typical main se-quence evolution. The stellar luminosity does not changeconsiderably from M ∗ = 30 M (cid:12) to 40 M (cid:12) - in fact there isa slight decrease. Contrary, for stars with mass >
40 M (cid:12) the stellar luminosity increases with mass. This is re-flected as a slight dip in the variation of Γ e with massbetween 30 M (cid:12) and 45 M (cid:12) . Thus, the values of Γ e de-rived for these two masses do not show considerable dif-ference (0.2381 and 0.2269 respectively). We concludethat the radiation force and, hence, also the dynamicsof the outflow do not change to a great extent for thesemasses.For simplicity we keep the inner jet launching radiusand the density at the base of the flow fixed for all mass Figure 11.
Same as figure 10, but for simulation runs for differentstellar mass. runs. The radiation force is not only altered by changesin Γ e , but also due to the physical scaling that appearsin the prescription of the force multiplier. In particular,this is the Keplerian velocity at the inner launching ra-dius l , which is naturally different for different massesand scales with √ M ∗ . One would then expect that theradiative force changes by a factor equivalent to √ M ∗ .The effect of such a change with increasing central massnot only enhances the mass flux launched from the diskboundary but also imparts the resulting outflow a widermorphology (see Table 2).In summary, with the above parameter study of differ-ent central stellar masses, we observe that as the lumi-nosity (or mass) of the central star increases, the stellarradiation force becomes relatively stronger. This leads toa higher degree of outflow de-collimation confirming theabove mentioned observational picture of outflow evolu-tion in massive young stars. Radiation field and jet mass flux
The simulations presented so far have been performedwith ”floating” boundary conditions at the base of theflow (see § ab initio , the disadvantage is thatthe mass flux emerging from the underlying disk bound-ary cannot be prescribed a priori . In fact, the mass fluxis self-consistently calculated each time step by ensuringcontinuity of outgoing waves between the domain andthe ghost cells. Thus, in our approach applied so far wefix the initial flow density and float the vertical outflowvelocity at the boundary.The total mass outflow rate ( ˙ M rad + ˙ M vert ) in physicalunits is shown in Fig. 12 for runs with different stellarmass and applying a floating boundary condition. Thetotal mass flux has a √ M ∗ dependence for pure MHDruns on converting it to physical units. While for asteady-state radiative MHD flow, the physical mass fluxalso depends on the Eddington parameter Γ e , which is6 Figure 12.
Outflow mass flux for simulations applying a differentstellar mass. The top panel shows the total mass outflow rates foreach star in case of a pure MHD flow ( stars ), and for the same out-flows including radiative forces ( dots ). The percentage differencebetween these mass fluxes is shown in the bottom panel . related to the stellar luminosity.The percentage change between these mass outflowrates is shown in the bottom panel of Fig. 12. Inter-estingly, the profile of this curves is similar to the varia-tion of the opening angles (see Fig. 11), indicating thatthe increased mass flux could in fact play a role for de-collimation. Thus, the combination of both the ”float-ing” boundary conditions and the disturbance of the gasphysics at the base of the flow by radiative source termsleads to the enhanced mass flux, which modifies the col-limation of the flow.However, also the direct effect of radiative force on theflow can physically deflect the flow and eventually lead tode-collimation. Thus, the problem becomes quite com-plex as the direct influence of radiative force on the flowis mingled with its second order effects for e.g. increasein the mass flux. In order to single out the influence ofdirect de-collimation effects by radiative forces, we haveperformed simulations where we have fixed the outflowmass flux as a boundary condition.It is essential to inject the outflow with a velocitywhich is supersonic to begin with. We fix the mass fluxat the boundary by setting the inflow velocity to v z =0 . × v Kep , resulting in a flow with slow magneto-sonicMach number close to unity and indicating a slightly su-personic flow. In Fig. 13 we show the mass outflow flowrates derived along the top ( z = z m ) and right ( r = r m )boundaries (in red and blue, respectively). We have cal-culated the radial and vertical mass outflow rates usingfollowing expressions,˙ M rad = r m (cid:90) z m ρv r | r m dz, (23)˙ M vert = (cid:90) r m rρv z | z m dr. (24)In the simulations with floating boundary conditions, Figure 13.
The variation of the outflow mass flux with time insimulation runs M30a055b3inj ( top panel ) and M30a055b5 ( bottompanel ). The red and blue lines show the variation of the verticaland radial mass fluxes, respectively (see Eqs. 23,24). The black line corresponds to the mass flux injected from the underlying disk. The green line is the total mass flux leaving the computational domain.Coinciding green and black lines proves conservation of the massflux. we observe that the mass flow rate reaches a steady-state,first derived self-consistently in the pure MHD flow, andthen suddenly increased on adding the radiative sourceterm. Essentially, this rise is not seen when we applyfixed-mass-flux boundary conditions. However, the factthat the curves of radial and vertical mass flux intersectafter ∼
60 disk rotations after switching-on the radiationforce, is a clear signature of de-collimation. Thus, themodification in terms of collimation due to the direct impact of radiative forces on a flow injected with constantmass flux can be now understood with our simulationruns with a priori fixed-mass-flux boundary conditions.We have thus performed three simulation runs with dif-ferent stellar mass applying the above-mentioned mass-injection boundary condition (see Table 2). The mostevident measure of de-collimation that can be observedis the ratio of vertical to radial mass flow rate. A higherratio indicates a higher degree of jet collimation. Thismass flux ratio for the chosen stellar masses of 25, 30,and 50 M (cid:12) changes from 1.12 after the steady-state pureMHD flow to 0.92, 0.85, and 0.80 respectively, when ra-diative force is considered. Note that the radiative forcesalso accelerate the flow, increasing the maximum veloc-ity for M30a055b3inj in steady-state radiative MHD bya factor of ∼ . direct influence of the stellar line-driven force on the out-flow, as the amount of mass flux launched into the do-main is fixed for these runs (it is thus not a 2nd-ordereffect resulting from a higher mass flux due triggered bythe radiation force).In summary, we conclude here that the line-drivingradiation force from the star plays a significant role inHD jets and radiation pressure 17directly modifying the dynamics of previously launchedMHD flow. This force not only aids in acceleration ofthe flow, but also pushes the outflow material away fromthe axis resulting in a substantially wider opening angle. Acceleration by the disk radiation field
We have also carried out a number of simulations thatinvolve only radiation forces caused by a (hot) accretiondisk. In principle, in addition to the intrinsic accretionluminosity, the disk radiative force can be enhanced byconsidering irradiation from the central star (Proga et al.1998, Drew et al. 1998). For the present model, however,we do not consider irradiation.Here we present results of two of our simulations, de-noted as Disk1 and Disk2 (see Table 4). We have ap-plied the parameter set as for the reference simulationM30a055b5, however, two more parameters governingthe dimensionless accretion rate that appears in the ac-cretion luminosity (Eq. B5) have to be prescribed. Thatis the ratio µ of disk to stellar luminosity, and the ratio Λof stellar radius to inner launching radius (see Table 1).Naturally, a disk extending to radii closer to the centerof gravity would have much higher temperatures thus lu-minosity (see Vaidya et al. 2009 for an application tomassive young stars) and the resulting radiation forceswould be able to affect the outflow more.For the number values applied for these additionalparameters we again follow the bloating star model(Hosokawa & Omukai 2009) assuming an accretion rateof 10 − M (cid:12) yr − , and an inner jet launching radius of l = 0 . µ = 0 . . ρ = 5 × − g cm − . As the radiation force isvery sensitive to the density, the lower density at the jetbase does increase the radiation force from the disk bya factor of three. Note that, physically, this increase ofthe disk radiation force could mimic the effects of stellarirradiation.We quantify the change in collimation degree in runsDisk1 and Disk2 again with the parameter ∆ ζ [%] (seeEq. 20). Simulation run Disk1 has a maximum of∆ ζ [%] = − .
8, while run Disk2 has ∆ ζ [%] = 3 . disk radiation force alone is not strong enough to de-collimate the flow .Interestingly, simulation run Disk2 shows some un-steady behavior close to the axis - a feature which isabsent when only stellar radiation forces were consid-ered. This is consistent with the findings of Proga et al.(1999), who detect a similarly unsteady behavior in theirsimulations for the case when the radiative force is dom-inated by the underlying disk. The fact that we do notsee this feature in simulation Disk1 suggests that the ra-diative force in run Disk1 is not comparable to the otherforces that control the flow dynamics, and that this out-flow does more correspond to a pure MHD flow. Theunsteady flow behavior in run Disk2 is also reflected inthe poloidal velocity evolution as steadily propagating ’wiggles’ in the velocity profile along the field line (orstreamline) v p ( s ) (Figure 6).As maximum outflow velocity for the simulations Disk1and Disk2 we obtain relatively high velocities of ∼ − . This is mainly due to the above mentioned factthat the outflow is launched deeper in the potential welland as close as l = 0.1 AU from the central star. Com-pared to the pure MHD run, we see, however, that the jetwhich is affected by disk radiation forces (only) achievesslightly higher asymptotic velocities, as the MHD reachesonly ∼
300 km s − . Thus, the disk radiation force affectsthe outflow primarily by slightly accelerating it. We donot see indication that the outflow collimation is affected,which is understandable since the disk radiation forceacts mainly in vertical direction.In summary, jet acceleration and collimation is ratherweakly affected by the disk radiative forces as their mag-nitude is orders of magnitude smaller compared to ra-diation forces from the central star. However, for otherastrophysical jet sources such as CVs and AGNs, the ra-diative force from the disk could play a significant role.Applying the scaling relations and comparing the amountof energy radiated per unit area from the standard thindisk from a typical young massive star (MS) and a stan-dard white dwarf (WD) we obtain from equation B1, D MS D W D = 0 . (cid:18) − M (cid:12) yr − g s − (cid:19) (cid:18)
30 M (cid:12) (cid:12) (cid:19) (cid:18) . cm (cid:19) − (25). This indicates that the efficiency of line driving due todisk force alone for such massive proto-stellar disks is lessthan the compact object accretion disks like CVs whichhave been the focus of previous studies by Pereyra &Kallman (2003). This ratio is even smaller if we considerdisks around AGNs as considered by Proga & Kallman(2004). Limitations of our model approach
The main goal of our study was to disentangle the ef-fects of radiative forces from the young star-disk systemon a pure MHD outflow launched within the standardpicture of magneto-centrifugal, magnetohydrodynamicjet formation. In this section we discuss the limitationsof our model approach to the subject of jet formationfrom massive young stars.
Prescription of radiation force
The prescription of radiation force used for the presentmodel does not explicitly consider the radiation transfer.Instead it implements the radiative force due to lines as asource term in the momentum and the energy equations(see § Possible existence of a stellar wind
Our model does not consider a stellar wind from themassive young star. Observed outflows around youngmassive stars are usually thought to have velocities ofabout 200 - 500 km s − , a value much smaller than thespeed of stellar winds in the evolved stages of OB stars (ofabout ∼ − ). This difference in velocity scalecould give a clue to difference in environment around thestar. During the formation stage, the star is surrounded8by dense gas and dust as compared to more evolved stagewhere it has cleared all surrounding matter. A rarerenvironment in a more evolved stage could lead to anefficient acceleration of winds via line driving to 1000km s − . Hence, for more massive and hotter young stars,impact due stellar winds could play a vital role on thedynamical evolution of disk winds. Lack of knowledge of system parameters
The inner regions around massive young stars are notaccesible with present day telescopes. Thus, many phys-ical parameters for these inner regions are not observa-tionally constraint (for e.g the inner disk radius, massdensity, magnetic field strength) In the present model, wefollow the notion that high-mass stars form as a scaled-up low-mass stars. Thus most of the parameters usedfor the present modeling are derived from estimates forlow-mass stars.
Ideal MHD and line force parameters
Simulations presented here are done using ideal MHDapproximation. This is in principle fine as ionizationsfractions are usually high enough to provide a good cou-pling between matter and field. However, for simplicitywe have assumed a constant ionization fraction through-out the wind. The radiative force parameters do have anexplicit dependence on the ionization fraction. The vary-ing temperature across the jet implies to varying degreeof ionization which would modify the line force parame-ters and thus the value of the force. In this model, we usean idealized assumption of a constant ionization fractionthroughout and fix the radiative force parameters for allsimulation runs. While in case of non-magnetic CVs,Pereyra & Kallman (2003) has shown that the qualita-tive effect of the line force is not altered by incorporat-ing local ionization effects for the line force parameters.However, the importance of such local ionization effectsfor the case of young massive stars is not known in theliterature.
Prescription of the disk dynamics
A self-consistent modeling would need to incorporatethe time-evolution of the disk structure in the simulation,and to treat accretion and ejection processes simultane-ously. Such simulations are currently performed in gen-eral applications of jet launching (e.g Casse & Keppens2002, Zanni et al. 2007, Murphy et al. 2010), however, itis too early to be applicable to massive young stars - onereason being the lack of knowledge of the accretion diskparameters (e.g. the question whether there is a thin orthick disk), another one that also radiative effects wouldhave to be implemented in these simulations. CONCLUSIONS
We have studied the impact of a line-driven radia-tion forces on the acceleration and collimation of a MHDjets,around young massive stellar object. Our main mo-tivation was to give a solid theoretical understanding forthe outflow evolution hypothesis presented by Beuther &Shepherd (2005). For the radiation forces we have con-sidered stellar and disk luminosity. Our basic approachwas to initially launch a MHD jet from the underlyingdisk surface, and then, after the pure MHD outflow has achieved a steady state, to switch on the radiative forcesand study their effect on the existing MHD outflow.We performed a number of simulations with the linedriving force exerted by stellar radiation. These simu-lations were performed for a wide range of physical pa-rameters, as i) the central stellar mass, ii) the magneticflux, and iii) the line-force parameter α . In order toapply our method of calculating the line driving force(CAK approach), we have modified the numerical codePLUTO to incorporate appropriate projections of gradi-ents of the 2 D velocity field along the light path. Addi-tionally, we have consistently implemented these projec-tions for different radiation sources properly consideringthe geometry of the radiation field. All these simulationshave ’floating’ and casually consistent inflow boundaryin which the mass flux is not fixed and is derived by thenumerical solution.Our main conclusions from this analysis are as follows.The line driven force from the central star for the pa-rameters considered does play a significant role in mod-ifying the dynamics in terms of collimation and acceler-ation of the outflow. We find that the outflow velocityis increased by a factor of 1.5 - 2 by radiation forces ascompared to the pure MHD flow. Also, the degree of col-limation is lowered, visible e.g. in a 30%-change in themagnetic flux profile, or the wider opening angle of themagnetic field lines.Investigating different stellar masses we determine theamount of de-collimation by measuring the opening an-gles of a typical field lines (that with the highest massflux) at the Aflv´ e n and fast magneto-sonic point. We findthat for a stellar mass of 20M (cid:12) , the opening angles are20 ◦ and 16 ◦ , respectively. For a 60M (cid:12) star these valuesincrease to 32 ◦ and 29 ◦ , indicating substantial amountof de-collimation due to the increased stellar luminosity.This findings confirm the observed evolutionary picturefor massive outflows in which more massive stars tendto have outflows of lower collimation degree. Note thatfor massive young stars, our results do not only consti-tute a relation of different stellar masses, but correspondalso to a time-evolution of outflow systems, as the centralmass can be substantially increased during massive starformation.We have also performed simulations with injection offixed mass flux from the disk boundary for various stellarmasses. We find that the ratio of vertical to radial massflux in these runs decreases from 0.92 to 0.80 with in-crease in stellar mass from 25 M (cid:12) to 50 M (cid:12) . This clearlyindicates the fact that the line driving force from cen-tral star plays an influential role in the physical under-standing of observed evolutionary picture pertaining tooutflows from young massive stars.So far, the magnetic field strength in the jet forma-tion region close-by a massive young star forming isan unknown quantity. We therefore have carried outsimulations with different field strength (plasma- β ∼ . , . , . α is verycritical in determining the magnitude of the line-drivenforces. Lower values of α leads to an efficient radiativeHD jets and radiation pressure 19force from the central star and thus decollimate the flowto a larger extent as compared to higher α values. Sincethe radiation forces also affect the mass outflow rates forsimulations, even small change in α may lead to signifi-cant changes in mass flux up to ∼ later stages of young massive star evolution. Inspite of these limitations, we find clear evidence of accel-eration and de-collimation of jets launched from massiveyoung stars.In summary, among all the radiative sources consid-ered to study the dynamics of outflow launched from theyoung massive star, we see that the line force from thecentral luminous star is very efficient in de-collimatingand accelerating the flow. The line force from the under-lying disk is not as significant as compared to the stellarforce. Also, dynamical effects on the outflow due to theoptically thin electron scattering continuum force fromthe central star and the disk is not significant. Further-more, we confirm the observed trend of de-collimationseen in outflows from massive stars at different evolu-tionary stages.We acknowledge the Klaus Tschira Stiftung for fund-ing this work carried out at Max Planck Institute ofAstronomy, Heidelberg and also convey our thanks tothe Heidelberg Graduate school of Fundamental Physics(HGSFP). We appreciate insightful comments fromDaniel Proga on this work and we thank him for his con-structive suggestions. We also like to thank A. Feldmeierfor his helpful comments for this work. We also forwardour thanks to T. Hosokawa for providing his data for thestellar evolution model considered here. APPENDIX
BASICS OF CASTOR, ABBOTT AND KLEIN THEORY
According to the Castor, Abbott and Klein theory (Castor et al. 1975), the force multiplier can be expressed as, M ( T ) = (cid:88) lines (cid:18) ∆ ν D F ν F − e − η T T (cid:19) , (A1)where ∆ ν D the Doppler shift, F ν the radiation flux at frequency ν , and F the total integrated flux. The optical depthparameter T is related to the gradient of velocity, the density in the wind and the ion thermal velocity v th , T = ρσ e v th | ˆ n · ∇ (ˆ n · (cid:126)v ) | . (A2)The optical depth parameter T can be related to the optical depth of a particular line τ L (ˆ n ) = η T . The line strength η is the ratio of the line opacity κ L to the electron scattering opacity σ e , while ˆ n is the unit vector along the line ofsight (l.o.s.). The classification of lines in optically thin and thick lines is done on the basis of interaction probabilities.Optically thick lines are those which have interaction probability of unity. Lines with optical depths τ L < τ L . Based on this approximation, the force multiplier is separableand can be written differently for optical depths very high or very low. When the gas is optically thick, τ L >
1, theforce multiplier depends only on the local dynamical quantities of the flow, M thick ( T ) = (cid:88) thicklines (cid:18) ∆ ν D F ν F T (cid:19) , (A3)while for the optically thin case τ L < M thin ( T ) = (cid:88) thinlines (cid:18) ∆ ν D F ν F η (cid:19) . (A4)In general, for a gas distribution with a mixture of optically thick and thin lines, the empirical form of the totalforce multiplier integrated over all lines can be expressed as a power law, M ( T ) ∼ k T − α , where k and α are line forceparameters.0 LINE DRIVING FORCE DUE TO DISK ALONE
In addition to the line radiative forces from the central star, we also take into account line forces from the underlyingdisk. The underlying disk cannot be considered as a point source but rather is an extended cylindrical source ofradiation. Further, the disk luminosity varies with the radial distance from the central star. In the present simulations,we consider that the underlying disk as a steady-state standard thin disk with a temperature profile given by Shakura& Sunyaev (1973). Thus, the energy radiated per unit area D ( r ) at cylindrical radius r is D ( r ) = 3 ˙ M acc GM ∗ πr (cid:34) − (cid:18) l r (cid:19) / (cid:35) , (B1)where ˙ M acc is the steady-state accretion rate onto a central star of mass M ∗ . The inner launching radius is l . In caseof line forces from the accretion disk, we calculate the radial and vertical radiation flux from the standard disk, S r and S z , similar to Pereyra et al. (2000). Further, the l.o.s. velocity gradient in the force multiplier applied in case of disksis for simplicity reduced to ∂v z /∂z , since the bulk of the radiation flux from the disk is in vertical direction, (cid:126)f line , disk = σ e c [ S r (cid:126)r + S z (cid:126)z ] M ( T ) . (B2)Here, S r and S z are the radial and vertical flux components. Both depend on the disk luminosity and are implementedin the code in their dimensionless form.In order to reduce to a dimensionless form, all the physical quantities should be written as a product of its value incode units with appropriate scale factor (see Eq. 12). S r and S z are radial and vertical components of the flux emittedfrom the disk surface.(Hachiya et al. 1998, Pereyra et al. 2000) - S r = (cid:90) π (cid:90) endl o D ( r (cid:48) ) π z α +1 ( r − r cos( φ )) (cid:104) ( r + r (cid:48) + z − rr (cid:48) cos( φ )) / (cid:105) α r (cid:48) dr (cid:48) dφ, (B3) S z = (cid:90) π (cid:90) endl o D ( r (cid:48) ) π z α +2 (cid:104) ( r + r (cid:48) + z − rr (cid:48) cos( φ )) / (cid:105) α r (cid:48) dr (cid:48) dφ, (B4)where D(r (cid:48) ) is the amount of energy radiated per unit area from the standard thin disk (Eq B1). In the presentsimulations, the accretion disk is treated as a boundary and the accretion of matter in the disk is not considered. Themass accretion rate that appears in the above flux radiative components of the disk has to be prescribed on basis ofdimensionless parameters. ˙ M acc = 4 πcl σ e Γ e Λ µ. (B5)The radiation flux given above has to be written in dimensionless form to be incorporated in the simulations. Theseflux components in the code units are given as follows - S r , code = 3 cGM ∗ πσ e l Γ e Λ µ (cid:90) π (cid:90) endl r (cid:48) c ) (cid:32) − (cid:115) r (cid:48) c (cid:33) z α +1c ( r c − r (cid:48) c cos( φ )) (cid:104) (( r c ) + ( r (cid:48) c ) + ( z c ) − r c r (cid:48) c cos( φ )) / (cid:105) α r (cid:48) c dr (cid:48) c dφ, (B6) S z , code = 3 cGM ∗ πσ e l Γ e Λ µ (cid:90) π (cid:90) endl r (cid:48) c ) (cid:32) − (cid:115) r (cid:48) c (cid:33) z α +2c (cid:104) (( r c ) + ( r (cid:48) c ) + ( z c ) − r c r (cid:48) c cos( φ )) / (cid:105) α r (cid:48) c dr (cid:48) c dφ. (B7)Using the above formulations, the dimensionless radial and vertical component of the line force from the disk canthen be obtained. Their respective contours are shown in Fig. B. f r , codeline , disk ( r c , z c ) = f rline , disk / ( GM ∗ /l ) = M c ( T ) S r , code , (B8) f z , codeline , disk ( r c , z c ) = f zline , disk / ( GM ∗ /l ) = M c ( T ) S z , code . (B9)HD jets and radiation pressure 21 Figure B1.