External photoevaporation of protoplanetary discs in sparse stellar groups: the impact of dust growth

1.6 MB·

by  Stefano Facchini

Mon.Not.R.Astron.Soc.000,1–18(2015) Printed29January2016 (MNLATEXstylefilev2.2)
External photoevaporation of protoplanetary discs in sparse stellar
groups: the impact of dust growth
Stefano Facchini1,2(cid:63), Cathie J. Clarke1 and Thomas G. Bisbas2,3,4
1InstituteofAstronomy,MadingleyRoad,CambridgeCB3OHA,UK
2Max-Planck-Institutfu¨rExtraterrestrischePhysik,Giessenbachstrasse1,85748Garching,Germany
3UniversityCollegeLondon,GowerPlace,London,WC1E6BT,UK
4DepartmentofAstronomy,UniversityofFlorida,Gainesville,FL32611,USA
6
1
0
Submissiondate
2

n
a ABSTRACT
J

7 Weestimatethemasslossratesofphotoevaporativewindslaunchedfromtheouteredge
2 of protoplanetary discs impinged by an ambient radiation field. We focus on mild/moderate
environments (the number of stars in the group/cluster is N (cid:38) 50), and explore disc sizes

] ranging between 20 and 250AU. We evaluate the steady-state structures of the photoevap-
R
orative winds by coupling temperature estimates obtained with a PDR code with 1D radial
S
hydrodynamicalequations.Wealsoconsidertheimpactofdustdraggingandgraingrowthon
.
h thefinalmasslossrates.Wefindthatthesewindsaremuchmoresignificantthanhavebeen
p appreciatedhithertowhengraingrowthisincludedinthemodelling:inparticular,massloss
- rates(cid:38) 10−8M /yrarepredictedevenformodestbackgroundfieldstrengths((cid:38) 30G )in
o (cid:12) 0
thecaseofdiscsthatextendtoR>150AU.Graingrowthsignificantlyaffectsthefinalmass
r
t lossratesbyreducingtheaveragecrosssectionatFUVwavelengths,andthusallowingamuch
s
a morevigorousflow.Theradialprofilesofobservablequantities(inparticularsurfacedensity,
[ temperatureand velocitypatterns) indicatethatthese windshavecharacteristic featuresthat

are now potentially observable with ALMA. In particular, such discs should have extended
1
gaseousemissionthatisdustdepletedintheouterregions,characterisedbyanon-Keplerian
v
rotationcurve,andwitharadiallyincreasingtemperaturegradient.
2
6
Key words: accretion, accretion discs — circumstellar matter — protoplanetary discs —
5
hydrodynamics—planetarysystems:formation
7
0
.
1
0
6 1 INTRODUCTION they are embedded. Even more interestingly, Mann et al. (2014)
1 (seealsoMann&Williams2010)haveshownthatthedustcom-
Thereisobservationalevidencethattheenvironmentofstarform-
: ponent of discs tends to be less massive in the immediate vicin-
v ingregionscansignificantlyaffecttheevolutionofprotostellarand
ity of the dominant O star in the ONC, θ1C. Mann et al. (2015)
i protoplanetarydiscs.Untilrecentyearstheonlyevidenceofsuch
X havenotobservedsuchatrendinanotherO-starsbearingcluster,
interplayhascomefromsingleobjectsshowingclearsignaturesof
r NGC2024.Formilderstarformingregions,westilllackobserva-
anon-goinginteractionbetweentheirdiscandtheambientenviron-
a tionalevidenceofhowimportanttheenvironmentcanbeinshap-
ment(e.g.O’Dell,Wen&Hu1993;Bally,O’Dell&McCaughrean
ing the structure and the evolution of discs. Upcoming spatially
2000; Vicente & Alves 2005). With the advent of sub-millimetre
resolvedsurveysofdiscsmightbeabletoclarifythis.Earlierstud-
surveys,wecannowstartaccessingstatisticalsamplesofrelevant
ieshaveanalysedwhetherdiscfrequencydependsondistancefrom
propertiesofprotoplanetarydiscs(mostlytheirmassandouterra-
themostmassivestarsinOBstars-bearingclusters.Mostofthem
dius)andinferwhethertheydoshowimprintsofbeingaffectedby
agreethatinthecloseproximitytotheluminousOstars,discfrac-
the environment. Such studies have shown that this is indeed the
tiondeclinesbyafactorof∼ 2.Thistrendhasbeenobservedfor
case in the most extreme star forming region within 500pc from
exampleinNGC2244(Balogetal.2007),inNGC6611(Guarcello
us,theOrionNebulaCluster(ONC).Inthiscluster,deJuanOve-
etal.2007,2009),inNGC6357(Fangetal.2012)andinCygOB2
laretal.(2012)havesuggestedthattypicaldiscsizesdecreaseasa
(Guarcelloetal.inprep.).However,theseresultsarestillcontro-
functionofthestellarsurfacedensityoftheenvironmentinwhich
versial:forexample,Richertetal.(2015)haverecentlyarguedthat
suchdecreasingtrendindiscfractioncouldbeenhancedbysample
incompletenessinsomeofthesestudies.Inlesspopulousclusters,
(cid:63) [email protected]
(cid:13)c 2015RAS
2 Facchini,ClarkeandBisbas
inwhichthenumberofsourcesnearOstarsissmaller,amorero- gravitational potential well of the central star. Since the outer re-
bustresultisthatnoevidenceofspatialgradientisobserved(e.g. gions of discs contain the bulk of the mass and since the rate of
Roccatagliataetal.2011). discevolutionisdeterminedbytheouterdiscradius,suchexternal
Therearetwomainenvironmentalmechanismsthatcanaffect winds have important implications for disc lifetimes and surface
protoplanetarydiscs:star-discinteractions,andphotoevaporation. densityprofiles.Indeed,asnotedbyClarke(2007)andconfirmed
Indenseclusters,duringthelifetimeofadisc(∼3−10Myr,e.g. byAnderson,Adams&Calvet(2013),suchwindsacceleratedisc
Haisch, Lada & Lada 2001; Mamajek 2009; Fedele et al. 2010), clearing not on account of the mass lost in the wind (which is a
gravitationalencounterswithotherstarsmayperturbthedisc,trun- fraction of that accreted onto the star over the disc lifetime) but
catingittoadefinedouterradius,andsteepeningthesurfaceden- insteadbecausetheymodifythedisc’souterboundary,preventing
sityprofile(seeHall1997;Breslauetal.2014;Rosottietal.2014; discspreadingandkeepingtheviscoustimescalerelativelyshort.
Vincke,Breslau&Pfalzner2015,andreferencestherein).Veryfew Lada & Lada (2003), Porras et al. (2003) and others have
objects that might be undergoing such encounters have been de- shown that the probability density function for cluster member-
tected(Salyketal.2014;Cabritetal.2006;Daietal.2015),since ship number N for clusters in the Solar Neighbourhood scales
thetimescaleofobservabilityisveryshort. with 1/N2. Therefore the cumulative probability for a star to be
Anarguablymoreimportantmechanismisthephotoevapora- born in a cluster of size (cid:54) N is inversely proportional to N
tioncausedbytheenergeticradiationpermeatingtheyoungasso- (e.g.Adams2010).Within2kpcoftheSun,themedianvalueof
ciations. When stars form in groups, EUV (Extreme Ultraviolet) thisdistribution(alsotakingintoaccountisolatedstars)is∼ 300
andFUV(FarUltraviolet)radiationofthemostmassivestarsheats (Adams et al. 2006; Allen et al. 2007). The majority of proto-
theouterregionsofprotoplanetarydiscs,andcandriveagaseous planetary discs are then very likely to evolve embedded in rel-
flow from the disc surface (Hollenbach et al. 1994; Adams et al. atively ‘mild’ (low N) environments. Clusters and groups with
2004, hereafter A04). Such a scenario has been very successful N (cid:54)500havetypicalFUVfieldsG (cid:54)3000G (e.g.Fatuzzo
FUV 0
in explaining the so-called ‘proplyds’ (Johnstone, Hollenbach & &Adams2008),whereG isthelocalFUVinterstellarfield(1.6×
0
Bally1998;Sto¨rzer&Hollenbach1999;Richling&Yorke2000), 10−3ergs−1cm−2,Habing1968).Fatuzzo&Adams(2008)have
i.e.darksilhouettediscsorcocoon-likestructuresobservedinstar alsoshownthattheFUVfluxdependsstronglyonthemostmassive
formingregions.Suchobjectshavebeenobservedforexamplein starinthecluster,andontheintra-clustercolumndensityofabsorb-
the ONC (O’Dell, Wen & Hu 1993; Bally, O’Dell & McCaugh- ingmaterial.Starsinsmallgroups(N ∼ 50−100)arelikelyto
rean 2000), in Cyg OB2 (Wright et al. 2012), in Carina (Smith, beimpingedbyanenvironmentalFUVfieldof∼30−300G ,i.e.
0
Bally & Morse2003), and in otherstar forming regions showing slightlylargerthanthelocalfield.Themajorityofprotoplanetary
massiveOBassociations(i.e.luminoussourcesofhighenergyra- discsarethenverylikelytolieinthesubcriticalregimeofexternal
diation).WhentheUVfluximpingingontothediscsislesssevere photoevaporation,atleastintheSolarNeighbourhood.
thaninthesestarformingregions,theeffectofsuchmasslosson Note that at such low values the external radiation becomes
theevolutionofprotoplanetarydiscs,andtheimpactontheirplanet comparabletotheradiationfromthecentralstarevenatlargeradii
formationpotential,ishoweverrarelyconsidered. (e.g.Berginetal.2003).Inparticular,Franceetal.(2014)looked
To be more quantitative, external photoevaporation can be atasampleofCTTSsandestimatedtheFUVfieldclosetothecen-
summarisedinthreedifferentregimes.WhentheEUVcomponent tralobjectbyderivingtheFUVcontinuumandhotgaslinespro-
isdominant,duetotheproximityofamassiveOstar(e.g.θ1Cin filesfromhighspectralresolutionHSTobservationsatFUVwave-
theONC),theEUVfluxreachesthesurfaceofthedisc,heatingit lengths.TheyobtainedFUVfieldsof∼107G at∼1AU,which
0
upto∼ 104 K.Acompletelyionisedflowisformed,drivinggas areduetobothstellarandaccretioncontributions.Byconsidering
outwardsandphotoevaporatingthedisc(Hollenbachetal.1994). geometrical dilution only, FUV fields of 300 and 30G are thus
0
When the EUV flux is less severe, the FUV field can generate a expectedat180and580AU,respectively.However,thebulkofthe
neutralwindthatisopticallythicktotheEUVradiation,sincethe discisshieldedbyitsinnerregions,especiallyforlowflaringan-
FUVcomponenthasalargerpenetratingdepththantheEUVone. gles,andmostoftheradiationisabsorbed∼2scale-heightsabove
Whenthisisthecase,twodifferentregimesmayoccur,depending themidplane,asshowninthermo-chemicalmodels(forexample,
onthediscsizeRd.WedefinethegravitationalradiusRg ofthe Bruderer et al. 2012, compute an FUV intensity of ∼ 100G0 at
discastheradiusatwhichthethermalenergyisequaltothegrav- ∼100AUandz/R∼0.2forHD100546).Besides,themidplane
itational binding energy (equivalently, as the radius at which the oftheouterdisccanbedirectlyimpingedbytheambientradiation,
soundspeedisequaltotheescapevelocity): andgascanberemovedmoreeasily.Thusintheouterregionsof
GM µm (cid:18) T (cid:19)−1(cid:18)M (cid:19) discsexternalphotoevaporationcouldplayadominant(oratleast
R = ∗ H ≈140AU ∗ , (1) comparable) role with respect to the photoevaporation driven by
g k T 1000K M
B (cid:12) the FUV radiation of the central star (Gorti & Hollenbach 2009;
whereT isthetemperature,M isthemassofthecentralstar,and Gorti,Hollenbach&Dullemond2015).Evidentlythisisanissue
∗
µthemeanmolecularweight(inthiscaseweusedatomicgaswith thatwillonlybefinallysettledusing2Dradiationhydrodynamical
µ = 1.3).WhenR > R ,thewindislaunchedsupersonically simulationsincludingbothinternalandexternalFUVsources.
d g
(e.g.Johnstone,Hollenbach&Bally1998).Suchdiscsaredefined The subcritical regime of external photoevaporation is the
tobesupercritical.However,moderateFUVfieldsarenotableto most difficult to probe observationally. However, there are some
heatuptheouterregionsofdiscstohighenoughtemperaturessuch observed signatures that might be indicating that such a mecha-
that R < R . A04 have shown that when this is the case, the nismisoccurringeveninquiteisolatedsystems(i.e.verylowex-
g d
flowstructurecanbedescribedbyanonisothermalParkerwind, ternalUVfluxes).ApotentialexampleofthisisthediscaroundIM
with the addition of a centrifugal term. This model uses the rea- Lup.Panic´etal.(2009)haveshownthatoutsidethemm-brightdisc
sonable assumption that the mass loss rate from the outer rim of (R =400AU;Pinteetal.2008),thegasstructureasprobed
out,dust
thediscdominatesthemasslossratefromthesurfaceofthedisc byCOemissionlinesindicatesasteepdecreaseinthesurfaceden-
(seetheAppendixinA04),wherethegasismoreembeddedinthe sityprofile,whichresemblestheprofilesobtainedbyA04intheir
(cid:13)c 2015RAS,MNRAS000,1–18
Externalphotoevaporationofprotoplanetarydiscs 3
summarise the photoevaporative model by A04, and explain how
ourdescriptiondiffersfromtheirs.InSections3-6wedetailthe
mainingredientsofourmodel;respectivelythethermalproperties,
gashydrodynamics,dusthydrodynamics,andtheiterationproce-
duretoobtainthefinalsolutions.InSection7wepresentourre-
sults,whicharethendiscussedinSection8.InSection9wesum-
mariseourconclusions.
2 COMPARISONWITHPREVIOUSWORK
Weconsiderthesubcriticalregime(R < R )previouslyinvesti-
d g
gatedbyA04andassumethatthediscisirradiatedbyanisotropic
FUVambientradiation.A04haveshownthattheratiobetweenthe
Figure1.Schematicillustratingaphotoevaporatingdiscimpingedbyan
masslossratefromthediscsurfaceM˙ andthemasslossfrom
caamnbbieendtiFreUcVtiornaadliaotriomno,rineitshoetrsoupbiccr,itdiceaplernedgiinmgeo(nRthde<soRurgce).oTfhierrraaddiiaattiioonn. the disc edge M˙ scales as M˙sur/M˙ ∼sur(Rd/Rg)1/2. Since they
Fromtheouteredgeofthedisc,asubsonicradialflowdevelops,untilit focus on the subcritical regime, where Rd < Rg, the mass loss
reachesthecriticalradiusRc,wheretheflowvelocityisclosetothesound rateisdominatedbytheradialflowemergingfromtheouterrimof
speed.Fromthispointoutwardsthevelocityofthewindisapproximately thedisc,andnotbythematerialflowingfromitssurface.Thusthe
uniform.Inthissubcriticalregime,themasslossfromthedisc’souterrim subsonicwindcanbedescribedwitha1Dradialmodel(seeFig.1
ismoresignificantthanthemasslossfromthediscsurface. foraschematicofthemodel).A04computedthetemperatureasa
functionoflocalgasdensitynandFUV(fromtheambientradia-
tion)opticaldepthτ byusingthephoto-dissociationregion(PDR)
flowmodels.O¨bergetal.(2015)havealsodetectedaDCO+double codebyKaufmanetal.(1999).Theythencoupledthistemperature
ringedstructure,whichcouldbetracingaradiallyincreasingtem- dependencewiththesteadystatemomentum/continuityequations
peraturegradientintheouterregionsofthediscinagreementwith inordertoiterativelyfindaself-consistentsteadystatesolutionof
thepredictionsfromexternalphotoevaporativemodels(thoughthe thegaseousflow.Themodelwepresentinthispaperisverysimi-
authorsinterpretthedataintermsofnon-thermaldesorptionofCO lartotheoneproposedbyA04,butcontainssomekeydifferences,
ices).Notethatsimilarfeaturesarealsopredictedbyphotoevapo- that will be described in more detail between Sections 3 - 5. In
rationmodelswherethewindsaredrivenbytheFUVfieldemitted particular:
bythecentralstarandbytheaccretingmaterial(Gorti&Hollen-
bach 2009; Gorti, Hollenbach & Dullemond 2015). From simple (i) Wetakedeviationsfromisothermalityintoaccountinlocat-
calculations of external photoevaporation models the gas flow in ingthecriticalpointoftheflow,incontrasttoA04whoimposethe
theouterregionsofthisdiscwouldbehighlydustdepleted,since conditionthattheflowistransonicatthelocationcorrespondingto
thedragforceisveryweakforthelargegrainsprobedbysubmm thesonicpointforisothermalgas.Thisself-consistentlocationof
observations, and such an effect is compatible with the dust de- thecriticalpointresultsinadifferentlocationandlocalflowveloc-
pletedouterregionsofthedisc.Thereareothersystemswherethis ityatthispointcomparedwiththeisothermalsolution.Ourresults
secondeffect(aradiallydecreasingdust-to-gasratio,asprobedby demonstratethatthisdifferenceisthemostimportantone,sinceit
line/continuumsizediscrepancies)isobserved(e.g.Pie´tu,Dutrey modifies mass loss rates significantly. Having located the critical
& Guilloteau 2007; Andrews et al. 2012; Rosenfeld et al. 2013; pointoftheflowweareabletointegrateinwardstothediscedge
deGregorio-Monsalvoetal.2013).NotehoweverthatBirnstiel& (in contrast to A04 who adopt an iterative scheme in integrating
Andrews(2014)haveshownthatsuchdiscrepanciescanalsobeex- outwardsfromthediscedge).
plainedbythesize-dependentradialdriftofdustparticles,without (ii) Theself-consistentlocationofthecriticalpointallowsusto
appealingtoanouterdiscwind.Wewillsuggestalternativeobser- obtainsolutionsoveralargerrangeofparameterspacethanA04.
vationaldiagnosticsthatcandiscriminatebetweenthetwoscenar- Specifically we are able to find solutions for lower values of the
ios. interstellarfieldGFUV (i.e.downto30G0)andforawiderrange
In this paper, we firstly aim to compute the mass loss rates ofouterdiscradii(outtoRd =250AU).
of discs affected by external photoevaporation in the subcritical (iii) Havingdeterminedaflowsolutionatfixeddusttogasratio
regime, byextending the parameter space investigatedby A04 to wethentakeaccountofthefactthatonlythesmallgrainsareen-
larger (but still subcritical) discs, and milder ambient fields. Our trainedintheflowandre-computetheflowsolutionimplementing
methodofsolutionbroadlyfollowsthatofA04butwithsomeim- thereduceddusttogasratiointheflow.
portantdifferencesrelatingtotheeffectofnon-isothermalityincor- (iv) Our code computing the temperature structure is different
rectlylocatingthecriticalpointoftheflow.Wealsoiteratetowards fromA04(seeSection3below).
aself-consistentsolutionwhichtakesintoaccountthefactthatonly
smallgrainsareentrainedintheflow.Notethattheeffectofpartial
entrainment of dust grains upon the thermal structure of the flow
3 TEMPERATUREESTIMATES-PDR
hasbeenrecentlyappliedtotheinternalphotoevaporationscenario
byGorti,Hollenbach&Dullemond(2015).Wefinallypresenttyp- FUVradiation(6 (cid:54) hν < 13.6eV)playsacrucialroleindeter-
icalradialprofilesofthemainhydrodynamicalquantities,andpro- mining the thermal structure and corresponding chemistry in so-
pose potential observational signatures of ongoing photoevapora- called photodissociation regions (PDRs) where gas undergoes a
tioninthissubcriticalregime. transition between ionised and molecular state. If a disc is inter-
The paper is structured as follows. In Section 2 we briefly acting with an ionising radiation field impinging externally, then
(cid:13)c 2015RAS,MNRAS000,1–18
4 Facchini,ClarkeandBisbas
10-18 10-18
3cm)1100--2109 tPoEt 3cm)1100--2109
erg/s/1100--2221 CHH_22i__ofpnohr_mdis erg/s/1100--2221 tCCo tIII
( (
e10-23 FUV_p e10-23 O I
rat10-24 ζCR rat10-24 CO
eating1100--2265 tcguharesb_mgr ooling1100--2265 gas_gr
H C
10-27 10-27
10-4 10-3 10-2 10-1 100 101 102 10-4 10-3 10-2 10-1 100 101 102
A A
V V
Figure2.Heatingandcoolingratesperunitvolumeforaslabofmaterialofdensityn = 104cm−3 impingedbyaFUVfieldofGFUV = 300G0,
afterithasreachedthermo-chemicalequilibriumcomputedwiththereducedchemicalnetwork.Thedust-to-gasratiois10−2andthePAH-to-dustratiois
2.6×10−2.Theheatingfunctions(labelledinthelegend)arerespectively:totalheatingrate(tot);photoelectricheating(PE);Cionisation(Cion);H2
formation(H2form);H2 photodissociation(H2phdis);FUVpumping(FUVp);cosmicrays(ζCR);turbulentheating(turb);chemicalheating(chem);
gas-graincollisions(gasgr).MoredetailscanbefoundinBisbasetal.(2012).
betweentheionisationfront(ifpresent)anditsdensegas,there-
Table1.Abundancesofspeciesusedinthepresentwork,fromAsplund
gionisopticallythintoFUVradiationbutopticallythicktoEUV etal.(2009).
radiation(hν (cid:62)13.6eV).Thiscausesthegastobeinastateofpar-
tialdissociationwherethetemperatureissetbytheFUVradiation H 4×10−1 Mg+ 3.981×10−5
field. Under these conditions, determination of the gas tempera- H2 3×10−1 C+ 2.69×10−4
turerequiresmodellingofvariousheatingandcoolingprocessesin He 8.5×10−2 O 4.898×10−4
termsofacomplexandnon-linearsetofdifferentialequationsde-
scribinganetworkofchemicalreactions.Inthelasttwodecades,
many groups have managed to develop numerical codes that are of 30, 300 and 3000G . The spatial extent of each density pro-
0
abletosolvesuchproblems(seeRolligetal.2007,forapaperde- file is chosen so that the visual extinction, A , is in the range
V
tailinginter-comparisonbetweensuchcodes).Belowwedescribe 10−7 (cid:54) A (cid:54) 10. Density is sampled every 0.1 dex, and A
V V
the3D-PDRcode(Bisbasetal.2012)thatweusetoobtainthegas every0.05dex.Theresultsaretheninterpolatedwithcubic-spline
temperature as a function of column density, grain opacity, local functionstoamuchfinergrid.Thecosmicrayionisationrate,ζ ,
CR
gasdensityandFUVfield. is taken to be ζ = 5×10−17s−1. The treatment of PAHs is
CR
discussedinthenextSection.Weevolvethechemistryineachsim-
ulationfor10Myr.Wehavecheckedthatthetemperaturebalance
3.1 3D-PDR isreachedonshortertimescales,∼ 10kyr.Aposteriori,wehave
verifiedthattheflowtimescaleisalwayslongerthan10kyr,thus
The3D-PDRcode(Bisbasetal.2012)isathree-dimensionaltime-
temperaturebalancewithintheflowisareasonableassumption.
dependentcodetreatingphotodissociationregionsofarbitraryden-
sitydistribution.Usinganiterationschemeitsolvesthechemistry
ateverycellconsistingofthegaseousstructureuntilthermaland
3.2 DustandPAHs
chemicalequilibriumhasbeenreached.3D-PDRusesasetofvar-
iousheatingandcoolingfunctionsfullydescribedinBisbasetal. Dustplaysakeyroleinthedeterminationofthegastemperature.
(2012),andshowninFig.2foraspecificcase.Thecodehasbeen Forexample,forthedensityrangeaddressedhere,themainheat-
used in various one-dimensional (Bisbas et al. 2014; Bisbas, Pa- ingmechanismofthegasisusuallyphotoelectricheatingfromthe
padopoulos&Viti2015)andthree-dimensional(Offneretal.2013, atomiclayersofPAHs(e.g.Weingartner&Draine2001;Croxall
2014;Gachesetal.2015)applications.The1Dversionofthecode et al. 2012), but it can also affect the temperature through other
UCL-PDRhasbeenbenchmarkedwithmanyotherPDRcodes(Rol- mechanisms such as H2 formation and gas-grain collisions. Sec-
ligetal.2007). ondly,dustcanhaveasignificantimpactonthechemistry,which
In this paper, we have adopted the code modifications de- eventuallysetstheabundancesofthemaincoolants,e.g.bycon-
scribed in Bisbas et al. (2014). We use a reduced version of the trollingtheamountofices.Notethatallthesemechanismsdepend
UMIST 2012 network (McElroy et al. 2013) which contains 33 onthetotalsurfaceareaofdustgrains.Finally,dustsetstheattenu-
species (including e−) and 330 reactions. Table 1 shows the ini- ationfactoroftheFUVradiationviathepenetratingdepthτ,where
tialabundancesusedbythe 3D-PDR codeatthebeginningofthe wesetτ = 1.8AV = NHσFUV andNH isthehydrogencolumn
calculations(takenfromAsplundetal.2009).Usingthefullchem- density.Adetaileddescriptionofthislasteffectwillbedescribed
icalnetworkof215speciescausesdifferenceinthetemperatureby inSection5.
upto∼10−15%.Sinceotherunknownparameters(suchaspoly- Typicalinterstellardustcanbemodelledwithasimpledistri-
ciclic aromatic hydrocarbon abundances, see below) lead to even butiondn˜/ds∝s−q,wheren˜isthenumericaldensityofdustpar-
largeruncertainties,weoptedforthereducednetworktosavecom- ticles,andsisthegrainsize.Itiswellknownthatintheinterstellar
putationaltime.Weconsiderone-dimensionaluniformdensitypro- mediumqassumesatypicalvalueq∼3.5(historicallylabelledas
fileswith102 < n < 108cm−3 interactingwithradiationfields MRNdistribution,fromMathis,Rumpl&Nordsieck1977).More
(cid:13)c 2015RAS,MNRAS000,1–18
Externalphotoevaporationofprotoplanetarydiscs 5
140 140
n=102
120 120 n=103
n=104
100 100
n=105
) 80 ) 80 n=106
K K
( (
T 60 T 60
40 40
20 20
0 0
10-4 10-3 10-2 10-1 100 101 10-4 10-3 10-2 10-1 100 101
A A
V V
250 250
200 200
150 150
) )
K K
( (
T100 T100
50 50
0 0
10-4 10-3 10-2 10-1 100 101 10-4 10-3 10-2 10-1 100 101
A A
V V
800 450
700 400
350
600
300
500
) )250
K K
400
( (200
T T
300
150
200
100
100 50
0 0
10-4 10-3 10-2 10-1 100 101 10-4 10-3 10-2 10-1 100 101
A A
V V
Figure3.GastemperatureasafunctionofvisualextinctionAV forthreedifferentvaluesoftheambientFUVfield:30,300and3000G0,fromtoptobottom.
Theleftpanelsareassociatedtoagrainsizedistributionwithmaximumgrainsizesmax = 3.5µm,therightpanelstoadistributionwithsmax = 1mm.
Linecoloursindicatedifferentgasnumericaldensitiesn,rangingbetween102−106cm−3(thelegendisreportedintherighttoppanel).
recently, interferometric observations have suggested a shallower etal.2014,andreferencestherein),tosizes(cid:38)1mm.Graingrowth
distributionintheopticallythinregionsofprotoplanetarydiscsat canaltertheheatingandcoolingprocessesmentionedabove.
submm - mm wavelengths, where q ∼ 2.5−3 (e.g. Ricci et al. Theonlyobservationindirectlyconstrainingthegrainsizedis-
2010a,b,2012),anddustevolutionmodelshavesuggestedthesame tribution within the photoevaporative winds is by Sto¨rzer & Hol-
qualitativeresult(e.g.Birnstiel,Ormel&Dullemond2011).More- lenbach (1999), who estimated the cross section in supercritical
over, the maximum grain size does increase substantially in the winds at FUV wavelengths from the location of the ionisation
midplane of protoplanetary discs (see the recent review by Testi front around proplyds in Orion. Their best estimate is σ =
FUV
(cid:13)c 2015RAS,MNRAS000,1–18
6 Facchini,ClarkeandBisbas
8×10−22cm2,wheretheerrorbarsareassumedtobeverylarge, tive grain growth, in this paper we prefer to keep the PAH abun-
sincetheusedsampleisverysmall(∼10objects)andtheestimate dance fixed, since there is no clear observational indication that
is model dependent. Sto¨rzer & Hollenbach (1999) and later A04 PAHswouldfollowthesmallgrains.
noticedthatsuchcrosssectionis∼ 0.3thetypicalISMone,thus InordertocomputethephotoelectricheatingduetoPAHs,we
indicatingmoderategraingrowthofthedustentrainedinthewind. adoptthetreatmentbyBakes&Tielens(1994),withtheadditional
We compute the cross section at FUV wavelengths (λ = modifications suggested by Wolfire et al. (2003) (with the φ
PAH
0.1µm) for an ISM-like dust distribution (q = 3.5 and s = factordefinedintheirequation21equalto0.4).
max
0.25µm, see e.g. Draine 2011) using the code described in Wy-
att & Dent (2002) (see their section 4.1). We use the optical
constants from Li & Greenberg (1997), for iceless silicate grains
with a porosity of 0.3. The absorption coefficients of the grains
Q (λ,s) are computed using Mie theory (Bohren & Huffman
abs
1983),Rayleigh-Ganstheoryorgeometricopticsintheappropri-
ate limits (Laor & Draine 1993). Assuming a dust-to-gas ratio to
3.3 Results
10−2 weobtainacrosssectionσ ≈ 2.6×10−21cm2,asex-
FUV
pected.Wethendeterminethemaximumgrainsizeofadistribution Fig.3showstheresultsofthe3D-PDRcodeforFUVfieldsof30,
withq=3.5thathasanassociatedcrosssection0.3timessmaller 300and3000G forthetwograinsizedistributions.Wepresent
0
thantheISM-likeone.Thebestfitgivessmax = 3.5µm.Thisre- thegastemperatureasafunctionofvisualextinctionAV,forlog-
sultconfirmsthattheestimatesbySto¨rzer&Hollenbach(1999)are arithmically sampled values of gas density n. As expected, the
probingmoderategraingrowth. gastemperatureincreaseswithintensityoftheexternalradiation.
InthePDRcode(andthroughoutthewholepaper),wecon- At a given G , for the parameters we have chosen, tempera-
FUV
sidertwoMRNdistributions(q =3.5)withmaximumgrainsizes tureisusuallyadecreasingfunctionofvisualextinctionA .For
V
s = 3.5µm and s = 1mm. We fix the dust-to-gas ratio s = 3.5µm temperature does also increase with A in the
max max max V
to10−2.Thefirstdistributionischosenastogivethesamecross marginally optically thick regime, when n ∼ 103 − 104cm−3.
sectionusedbyA04intheirmodels,inordertohaveadirectcom- Similarly,temperatureisnotamonotonicallydecreasingfunction
parisonwiththeirresults.Theseconddistributionrepresentsadisc ofdensityn.Thisnon-monotonicbehaviourisawellknownresult
where substantial grain growth has occurred, and the population inPDRcodes(Kaufmanetal.1999).
ofsmallgrains,whichcontributethemosttothedustsurfacearea Theseresultsdependsonthemetallicityandonthedust-to-gas
whenq>3,isreducedbyafactor∼102. ratio,sincetheybothregulatetheheatingandcoolingfunctions.In
PolycyclicAromaticHydrocarbons(PAHs)arethedominant particular, reducing the amount of small grains has the net effect
heatingsource,sincephotonsareveryefficientinremovingelec- ofslightlydecreasingthetemperatures,aswecanseecomparing
tronsfromsuchmono-layeredmolecules.However,theamountof thetemperaturesobtainedwiththetwodifferentdustdistributions.
PAHs in protoplanetary discs is not well constrained by observa- The smaller dust total surface area reduces some heating mecha-
tions. Infra-Red observations from Spitzer and from the ground nisms,suchasphotoelectricheating,H formationrateandH pho-
2 2
haverevealedstrongPAHemissionfeaturesinsomeHerbigAe/Be todissociation.However,thetotalheatingisstilldominatedbythe
stars, whereas very few T Tauri stars have shown them at a de- photoelectriceffectonPAHs.Coolingisnotaffectedsignificantly,
tectablelevel(Geersetal.2006,andreferencestherein).Thisdif- sincethemaincoolantsareCIIandOIintheopticallythinregime
ference is probably due to the more intense UV field emitted by (seeFig.2).Theneteffectisthattemperatureisnotaffectedsignif-
theAe/Bestars(e.g.Visseretal.2007)excitingthePAHs.Forthe icantlybythetotalamountofsmallgrains,providedthatthePAH
samereason,PAHemissionisspatiallyconcentratedintheinnerre- abundanceiskeptfixed.Theonlyexceptionisapparentinthehigh
gionofthediscs,wheretheUVfluxfromthecentralstarishigher densityandhighFUVfluxcase(n=106andG =3000G ).
FUV 0
(Geersetal.2007;Maaskantetal.2014).PAHabundanceindiscs In this regime, heating from FUV pumping of H molecules be-
2
isstillbeingdebated.Observations(e.g.Geersetal.2006;Oliveira comescomparabletothephotoelectricheating(Hollenbach&Tie-
etal.2010)tendtoindicatethatPAH-to-dustmassratioindiscsis lens1999,A04).TheabundanceofH moleculesdependsonthe
2
∼10%ofthePAH-to-dustmassratiointhegalacticISM,whichis H formationrate,whichdependsonthetotaldustsurfacearea.By
2
≈4×10−2(Draine&Li2007;Tielens2008),buttheestimateis reducingtheamountofsmallgrains,H formationratedrops,and
2
poorlyconstrained.Inthermo-chemicaldiscmodels,differentval- sodoestheFUVpumping,resultinginalowertemperature.
uesofPAH-to-dustmassratiohavebeenusedintheliterature,e.g. For the cases considered here and displayed in Fig. 3, at
0.5×10−2 (Bruderer et al. 2012), or 1.3×10−2 (Woitke et al. fixed density the temperature always shows a plateau at low op-
2015). Because of the level of the uncertainties in this value, we tical depths, and then decreases steeply with visual extinction,
simplyfixthePAH-to-gasmassratioto2.6×10−4,asinWolfire until it approaches unity and the gas becomes optically thick to
etal.(2003),sinceweusetheirprescriptiontocomputethephoto- theexternalradiation.Thesetwocharacteristicswillshapethefi-
electricheating.Forthedust-to-gasratioof10−2usedinthispaper, nal profiles of the subsonic winds. Finally, the temperatures are
thisvalueyieldsaPAH-to-dustmassratioof2.6×10−2. lowerthantheonescomputedbyA04(e.g.byafactor∼ 2when
Another debated topic is whether grain growth affects the G = 300G ),intheregionsofparameterspacethatoverlap,
FUV 0
amount of PAHs, i.e. whether they should scale according to the eveninthehottestcasewhens = 3.5µm.Forthesamedust
max
amount of small grains. This is still largely unconstrained. As an distribution, we obtain temperatures similar to the ones by Kauf-
example,oneofthestrongestPAHemissionfeaturesisobservedin manetal.(1999)(andthereforebyA04,whousethesamePDR
OphIRS48(Geersetal.2007),wheregraingrowthhascertainly code)byusingtheprescriptionbyBakes&Tielens(1994)(with-
occurred (e.g. van der Marel et al. 2013). Even though some re- outthecorrectionbyWolfireetal.2003)forthePAHphotoelectric
centmodels(e.g.Gorti,Hollenbach&Dullemond2015)assumea heating. In particular, we reproduce the same trend observed by
reducedPAHabundancewhensmallgrainsaredepletedbyeffec- Kaufmanetal.(1999)intheirfig.1.
(cid:13)c 2015RAS,MNRAS000,1–18
Externalphotoevaporationofprotoplanetarydiscs 7
4 HYDRODYNAMICS ofthecentrifugalterm.Bysolvingthisequation,wecanobtainthe
velocitystructureoftheflow,fromwhichwecanobtainthedensity
In the last section we have shown how to derive the temperature
structureofthegasviathecontinuityequation.
ofthegasinthePDRregionasafunctionoflocaldensitynand
visualextinctionA .Wenowwanttocoupletheseresultstothe
V
hydrodynamicalequationsdescribingtheflowdepartingfromthe 4.2 Thecriticalradius
edgeofaprotoplanetarydisc.
Thestructureofequation7showsanaturaldefinitionofacritical
point:whenther.h.s.oftheequationequals0.Whenthishappens,
4.1 Equationsforthewindstructure thel.h.s.ofthesameequationhastoequal0aswell.Thiscanhap-
penwhenuhasanullgradient,orwhenthesecondmultiplicand
We develop a 1D description of the problem in a quasi spherical
is0.Werequirethissecondcasetobetheone,sincewearelook-
geometry, where we self-consistently solve for the radial steady-
ingfortheanalogueofthetransonicsolutioninthetypicalParker
statestructureofthewindlaunchedfromthediscouteredge.The
wind problem (Clarke & Carswell 2007). In the isothermal case,
equationsshowninthissectionarebasedonthesimilarmodelby
this same procedure defines the sonic radius R , where u2 = f.
A04. s
Weexpectthecriticalradiustobeofthesameorderofmagnitude
Thetimeindependentversionofthecontinuityequationis:
asthesonicradius.Byrequiringthatthel.h.s.oftheequationis0,
M˙ =4πR2FµmHnv, (2) butdu/dξ(cid:54)=0,weimplicitlydefineacriticalvelocity:
wherem isthehydrogenatommassandvisthevelocityofthe ∂f
H u2 =f +g , (8)
gas.F isthefractionofthesolidanglesubtendedbythediscouter c ∂g
edge (see equation 14 by A04). The flow might be not perfectly
whichisrelatedtothesoundspeedu viaanadditionaltermtaking
s
spherical, since the wind might not be pressure confined at the
intoaccountpossibledeparturesfromisothermality.Theequation
boundariesofthewedgedefinedbythesolidangle4πF.Wecould
forthedimensionlesscriticalradiusξ =R /R isthen:
c c d
modeltheflowwithahyper-sphericalcontinuityequation,where
M˙ ∝ Rαandα > 2(α = 2definesthesphericalcase,seeequa- gτ ∂fξ3+2u2ξ2−βξ +β =0, (9)
tion2).Forsimplicity,weconsiderthesphericalcaseonlyinthis d∂τ c c c c
paper. where all the quantities are evaluated at the critical radius. This
Similarly,thetimeindependentmomentumequationis: problemhastobesolvedimplicitly,sinceallthequantitiesdepend
dv 1dP GM j2 onwherethecriticalradiusislocated,andontheinitialconditions
v + + ∗ − =0, (3) oftheproblem.Locatingthecriticalradiusself-consistentlyisakey
dR ρdR R2 R3
ingredientofourapproachinordertobeabletointegratetheflow
where j2 = GM∗Rd. We are considering a flow with uniform structureinwardsfromthispointtothediscouteredge.Notethat
specific angular momentum, given by the Keplerian angular mo- this is a key difference between our model and the one proposed
mentumattheouterradiusofthedisc.Weassumeanidealgaslaw by A04. In Section 7 we will show that taking into account this
forthepressureP = nkBT,whereTisafunctionofdensityand non-isothermaltermsignificantlyaffectsthefinalmasslossrates.
opticaldepth(asshowninSection3). ForeverysetofparametersM ,R andG thereisafam-
∗ d FUV
FollowingA04,wecanexpresstheaboveequationsinadi- ily of solutions each of which has a different value of T , mass
c
mensionlessform,whereourparametrisationchoiceisslightlydif- lossrateandpressureatR .Werequirethatourflowisinpressure
d
ferent. We define ξ ≡ R/Rd, f ≡ T/Tc, g ≡ n/nc and u ≡ equilibrium with the disc at Rd; therefore in principle we could
v/cs,c,wherecs isthesoundspeedofthegas.Withthesubscript usethepressureatthediscedgeasafurtherindependentvariable
c we indicate quantities evaluated at a critical radius Rc > Rd, whichwoulddefinethemasslossrateandflowstructureforagiven
whichisdefinedinSection4.2.Notethatwiththesedimensionless discinagivenenvironment(cf.A04).Thishoweverrequiresiter-
unitsthelocalsoundspeedu2s iscoincidentwiththetemperature ativesolutionsforvariousguessedflowvelocitiesatRdandisnot
f.Wealsodefinetheparametersβ: apreferredmethodinthecasewherethecriticalpointconditions
GM µm GM arethemselvesafunctionoftheflowsolution.Sincewestarteach
β = ∗ H = ∗ , (4) flowsolutionfromtheself-consistentlydeterminedcriticalpointof
R k T R c2
d B c d s,c
theflow,itisoperationallyconvenienttouseT astheindependent
c
andtheopticaldepthτ ≡NHσFUV,where variable.Eachsolutioncanthenbemappedontothecorresponding
(cid:90) ∞ valueofthepressureatthediscouteredge.
N = n(R(cid:48))dR(cid:48) (5)
H InordertoevaluatethecriticalradiusforagivenvalueofT
c
R
weproceedasfollows.Wefirstguessthevalueofthecriticalradius
isthecolumndensitybetweenRandinfinity.Wethusobtain:
(i.e.asparametrisedbyξ );foraninitialguesswetakethere-scaled
c
dτ sonicradiusξ ,usingtheisothermallimitofequation9:
=−σ R n g=−τ g, (6) s
dξ FUV d c d
β(cid:20) (cid:18) 8(cid:19)1/2(cid:21)
where τd = σFUVRdnc. Combining the continuity and the mo- ξs = 4 1+ 1− β , (10)
mentumequation,weobtainasingledifferentialequationforthe
dimensionlessvelocityu: andnotingthatβisfixedforgivenTcbyequation4.
Once we have an initial guess for the critical radius, for ev-
dlnu ∂f 2 ∂f ξ−1 ∂f
dξ (u2−f −g∂g)= ξ(f +g∂g)−β ξ3 +τdg∂τ.(7) eryiterationonthevalueofξc wecalculatethepairofvaluesfor
the number density and optical depth at this point for which the
ThislastequationreducestothestandardParkerwindequa- PDRmodelspredictthatT =T locally.Thissteprequiresanas-
c
tion(Parker1958),inthelimitofisothermalityandintheabsence sumptionaboutthedensitystructureoftheflowoutsidethecritical
(cid:13)c 2015RAS,MNRAS000,1–18
8 Facchini,ClarkeandBisbas
450
3.8
400 30G0
300G
350 0 3.6
3000G
0
300 U)3.4
K)250 /A
( Rc3.2
Tc200 g(
150 lo3.0
100
2.8
50
0
10-4 10-3 10-2 10-1 100 101 20 40 60 80 100 120 140
τ T (K)
c c
Figure4.Thecriticaltemperature,i.e.theinputparametertoobtainaso-
3.8
lutiontothethermo-hydrodynamicalequations,asafunctionoftheop-
ticaldepthevaluatedatthecriticalradiusRc,forGFUV = 30,300and
3.6
3000G0 forsolutionwithRd = 120.Thisimplicitrelationdefinesthe
bdoasuhneddarlyinceosntdoitsiomnasxat=R1=mmR.cE.vSeorylidculirnveesirseffeorrtoonesmspaexcifi=cs3o.5luµtimon, AU)3.4
soonlluyt,isoinnsce(sietedSepeecntidosno5n.1t)h.evalueofσFUV,whichvariesfordifferentflow R/c3.2
(
g
o
l3.0
radius.HerewefollowJohnstone,Hollenbach&Bally(1998)and
A04inassumingthatthevelocityisapproximatelyconstantinthis 2.8
regionsothatfromthecontinuityequationwethenobtain:
n(R)=n (cid:18)Rc(cid:19)2, forR>R , (11) 20 40 60 80 100 120 140
c R c T (K)
c
andconsequently
(cid:90) ∞ (cid:18)R (cid:19)2 Figure5.Toppanel:criticalradiusdependenceoncriticaltemperaturefor
τ =σ n c dR(cid:48) =σ R n . (12)
c FUV c R(cid:48) FUV c c GFUV = 30G0andsmax = 1mm,forasetofdiscradiirangingfrom
Rc 20AU(topline)to250AU(bottomline),sampledevery10AU.Every
Figure 4 gives an example of Tc as a function of τc for an lineillustratesadifferentdiscradiusRd.Bottompanel:solidlinesarethe
initialguessRc ∼500AUandforarangeofFUVfluxes.Forlow sameasinthetoppanel,forRd = 20,100and180AU(blue,blackand
FUVfluxes,thearenosolutionsinthelowopticaldepthrange(τ < greenlines,respectively).Dashedlinesindicatethesonicradius,asdefined
10−3).Thisconstraintissetbythelowerlimitwehaveimplicitly byequation10.
imposedonthedensity,whenwehavenotdealtwithanythingless
densethan102cm−3inSection3.
Wearenowabletodetermineτ andn bysolvingtheequa-
c c
4.2.1 Extendingthesolutions
tion:
T −T(τ )=0. (13) Inthedefinitionofthesonicradius(seeequation10),byconstruc-
c c
tionwearerequiringthatβ >8.Thusweareimplicitlysettingan
We can evaluate the partial derivatives of T with respect to upperlimittotheT wecansetatagivenR :
c d
densityandopticaldepth(atfixedR )andusetheseinequations8
c
and9.WethensolveforthenextvalueofRc bysolvingequation T < GM∗µmH (cid:39)873(cid:18)M∗(cid:19)(cid:18)20AU(cid:19)K. (14)
9usingabisectormethod.Wetheniteratethewholeprocess,until c 8k R M R
B d (cid:12) d
we reach convergence on the value of ξ . By doing so, we have
c
obtainedthecriticalradius,andwehavetheinitialconditionsfor Thisstronglylimitstheexplorableregionofparameterspace,when
T (givenbytheinputparameterT ),v,nandτ.Thecriticalradius R islargeandtheFUVflux(andthereforetemperature)ishigh.
c d
isalwayslargerthanthesonicradius(asdefinedbyequation10), In order to enlarge the parameter space, when β < 8, as a first
butatthemostbyafactorofafewtenthsofadex(seeanexample guessforR wechoosethecriticalradiusofthesolutionwiththe
c
inFig.5).Ontheotherhand,thecriticalvelocityisusuallysmaller same T and the closest value of R showing a iteratively con-
c d
thanthesoundspeedatthecriticalradius,by∼ 20%atthemost vergedsolutionforthecriticalradius.InFig.5weshowR versus
c
(seeanexampleinFig.6).Finally,wecanestimatethemass-loss T when G = 30G , for a range of disc radii R between
c FUV 0 d
rateM˙ fromequation2.Wearereadytostartintegratingequation 20 and 250 AU, sampled every 10 AU. The critical radius does
7inwards,fromthecriticalradiustothediscouteredge,inorder notshowastrongdependenceonR ,atfixedT .Inthiswaywe
d c
todeducethevelocityprofileofthesubsonicflowbetweenR and areabletoobtainmoresolutions,eventhoughinsomeregionsof
d
R . parameterspacewedonotmanagetoobtainone.
c
(cid:13)c 2015RAS,MNRAS000,1–18
Externalphotoevaporationofprotoplanetarydiscs 9
bothradiusandvelocity:ξ=ξ +δξ,andu=u +δu.Weobtain
c c
1.05
thefollowingrelation:
√
(cid:12)
1.00 δu(cid:12)(cid:12) = −B+ B2−4AD, (15)
δξ(cid:12) 2A
ξc
0.95 where:
c
cs, 2g∂f g2∂2f
v/c0.90 A=2uc+ u ∂g + u ∂g2; (16)
c c
∂f 8g∂f 4g2∂2f ∂2f
0.85 B =2τdg∂τ + ξ ∂g + ξ ∂g2 +2τdg2∂τ∂g; (17)
0.80 (cid:18)2u2 2ξ−3 8g∂f 4τ g∂f
20 40 60 80 100 120 140 D=uc× ξ2c −β ξ4 + ξ2 ∂g + ξd ∂τ
T (K)
c 4g2∂2f 4τ g2 ∂2f ∂2f(cid:19)
+ + d +τ2g2 , (18)
ξ2 ∂g2 ξ ∂τ∂g d ∂τ2
Figure6.RatioofcriticalvelocitytothesoundspeedatRcforsolutions where all the quantities are evaluated at ξ = ξ . In equation 15
withGFUV=30G0andsmax=1mm,andRd=20,100and180AU c
wechosethepositiveroot,inordertopickthesubsonicsolutionat
(blue,blackandgreenlines,respectively).
ξ<ξ .
c
SimilarlytoMurray-Clay,Chiang&Murray(2009),wecal-
4.3 Methodofsolution culatethevelocitybyusing:
(cid:12) (cid:12)
Oncewehavelocatedthecriticalradiusanddeterminedthebound- du du(cid:12) δu(cid:12)
aryconditionforflowdensity,opticaldepthandvelocityforagiven dξ =Fexactdξ(cid:12)(cid:12) +(1−Fexact)δξ(cid:12)(cid:12) , (19)
T ,adiscradiusandaFUVfluxG ,wecanintegrateequation7 exact ξc
c FUV
fromRctoRd(inrescaledunits,fromξcto1).WeuseasimpleEu- wheredu/dξ|exactisevaluatedfromequation7,and
lercode.Ateverystepwecomputeui+1fromequation7.Then,we (cid:20) (cid:18) f g ∂f(cid:19)(cid:21)
calculatedimensionlessdensitygi+1fromthedimensionlessform Fexact =−erf h 1− u2 − u2 ∂g , (20)
ofequation2,andτ fromequation6.Atgiveng andτ
i+1 i+1 i+1
wecancomputethenewtemperaturepartialderivatives∂f/∂g where erf is the error function. The parameter h = 20 gives
i+1
and∂f/∂τi+1fromthePDRoutput,andthereforeobtainthenext the smoothing length of the transition between du/dξ|exact and
value of the temperature fi+1. With this set of equations we can δu/δξ|ξc whencomputingu(ξ).Weusethismodifiedversionof
self-consistentlysolvefortheflowsteady-state.Weuseaninitially equation7untilFexact =1tothelevelofmachineprecision.
logarithmicallyspacedgrid,sincetheabsolutevalueofthegradient
ofmostofthequantitiesincreasesasthesolutionapproachesR .
d 4.3.2 TemperaturecorrectionsnearR
Thishappensmostlyinthetransitionbetweentheopticallythinand d
theopticallythickregime,whenthetemperaturedropssteeply(as Somesolutionsbecomecompletelyopticallythickbeforereaching
showninFig.3).Inordertobetterresolvethisregionwealsoapply thediscouterradius(i.e.R > R ).InthePDRcalculations,we
d
anadaptivemeshalgorithm,i.e.werequirethatδξismallenough havesetaminimumtemperatureequalto10K.However,thetem-
toensurethattherelativechangeinvelocitybetweentwostepsis peratureoftheflowcouldbehigherthanthat,duetotheimpinging
lessthan1%(|ui+1−ui|/ui <0.01). radiationfromthecentralstar.Whenwecomputethetemperature
However, we need to apply slight modifications to the algo- oftheflow,wethereforeincludeheatingfromthecentralstar,by
rithmatthetwoboundariesoftheintegration,intheproximityof usingthefollowingsimpleprescription:
boththecriticalradiusR andthediscouterradiusR .
c d
T =max(T ,T ), (21)
PDR rad
whereT isthetemperatureevaluatedviathePDRcode(i.e.the
4.3.1 Thecriticalpoint PDR
oneusedsofarinthepaper),andT isgivenby:
rad
InSection4.2wehavedefinedthecriticalradiusastheradiusat
(cid:18) R (cid:19)−1/2
whichther.h.s.ofequation7vanishes.Moreover,atthissamera- T =100K , (22)
rad 1AU
dius the second multiplicand of the l.h.s. of the same equation is
0.SincewearestartingintegratingthesameequationfromRc,we i.e.atemperatureprofilefoundtofitthespectralenergydistribu-
have a null value both at the numerator and at the denominator. tionsofpassivediscs(e.g.Andrews&Williams2005,2007).
Thisleadstothepossibilityofhavingmultipletranscriticalsolu-
tions (the analogue version of the transonic solutions in the stan-
dard isothermal case). More specifically, we will have a solution
5 DUSTCOMPONENT
thatwillbesupersonicbetweenR andR ,andonesolutionthat
d c
willbesubsonic.Wearelookingforthesecondone.Thereforewe InSection4,wehavereportedtheequationsandtheprocedureto
needtoenforcethesolutiontorelaxontothesubsonicbranch.We obtainasolutionforthegasquantitiesbetweenthecriticalradius
dosobyfollowingaprocedurethatissimilartotheoneusedby and the disc’s outer edge. However, this solution depends on the
Murray-Clay,Chiang&Murray(2009)forananalogousproblem. dustpropertieswithintheflow.Inparticular,asmentionedinSec-
Weexpandequation7aroundthecriticalpointtofirstorderin tion3.2,theattenuationfactorstronglydependsonthegrainsize
(cid:13)c 2015RAS,MNRAS000,1–18
10 Facchini,ClarkeandBisbas
distributionandonthemaximumgrainsizes entrainedinthe
entr
wind.SuchdependencecanbesummarisedinhowσFUVisrelated 10-21
tothesetwopropertiesofthedustmaterial.Inthissection,wealso
describethehydrodynamicequationsforthedustparticles,which
willbeusedtodeterminethemaximumgrainsizeentrainedinthe )10-22
flow,andthusthecrosssection. 2m
c
(10-23
V
5.1 Crosssection U
F
σ
Itiswellknownthatprotoplanetarydiscswitnesssubstantialgrain 10-24
growthinthediscs’midplane(seeSection3.2).Suchgraingrowth
canbeschematisedintwoeffects:producingashallowerdistribu-
tion, i.e. a lower q, and leading to a larger maximum grain size 10-25
s .AsintroducedinSection3.2,inthispaperwefocusontwo 10-7 10-6 10-5 10-4 10-3 10-2 10-1
max
distributionswiththesamepowerlawindexq = 3.5butdifferent s (cm)
entr
maximumgrainsizes(s =3.5µmand1mm).Thesamedistri-
max
butionshavebeenusedtocomputethedependenceofthetempera-
tureonnandA .Asmentionedabove,weconsideradust-to-gas Figure7.Crosssectionatλ=0.1µmasafunctionofthemaximumgrain
V
ratioof0.01forboththedistributions. size entrained in the flow, for smax = 1mm (black line) and smax =
We compute the cross section at FUV wavelengths (λ = 3.5µm (red line). The initial distribution assumes a dust-to-gas ratio of
0.01. The cross section does not vary by much when the distribution is
0.1µm)usingthecodeofWyatt&Dent(2002)asbrieflydescribed
inSection3.2.Intheequationsreportedbelow,wedescribeagen- truncatedtosentr(cid:38)λ,sincethelargestcontributiontothecrosssectionin
thegeometriclimitcomesfromthesmallestgrainswhenq>3.
eral treatment where q can assume different values. We take into
accountthatthereisamaximumgrainsizeentrainedinthephoto-
evaporativewind,andwecomputethecrosssectionsforthetrun-
particles. We discretise the grain size distribution in N = 50
bin
cateddistribution:
bins,wherethegrainsizeissampledlogarithmicallybetweens
min
q
σFUV =µmHm sδ , (23) andsmax.Thedustmassdensitycanbewrittenas:
entr gd
4
where ρj (cid:39) n˜j3πρ¯s3max,j, (27)
(cid:90) sentr
q = πs2Q (λ,s)s−qds, (24) wherethesubscriptj referstothej-thbinofdustparticles.Here
s abs
smin ρj isthemassdensityofthej-thbin,andsmax,j isthemaximum
4 πρ¯ grainsizeofthej-thbin.
m = (s(4−q)−s(4−q)), (25)
entr 34−q entr min Thesteadyversionofthemomentumequationfordustparti-
clesis:
and
(cid:32) (cid:33) dv GM j2 F
s(4−q)−s(4−q) v j + ∗ − + D =0, (28)
δ =100 max min . (26) jdR R2 R3 m
gd s(4−q)−s(4−q) j
entr min
wherewehaveassumedthatdustispressureless,andthatthebuoy-
Intheseequations,δ istheeffectivegas-to-dustratiowithinthe
gd ancytermrelatedtothegaspressuregradientisnegligible(which
flow,andρ¯isthemeanmassdensityofadustgrain(inthispaper
istypicallythecaseforastrophysicalsystems,seee.g.Riceetal.
1g/cm3).Theminimumgrainsizeofthedistributionhasbeenset
2006;Laibe&Price2012).Thefluidquantitiesthatdonotshow
tos =5×10−7cminthewholepaper.
min anysubscriptrefertogasproperties,asinSection4.Thequantity
ThecrosssectionsarereportedinFig.7,asafunctionofthe
m is the mass of the j-th dust particle. The last term F repre-
j D
maximumgrainsizeentrainedintheflow,fors =3.5µm(red
max sentstheaerodynamicforceexperiencedbythedustparticle.This
line) and s = 1mm (black line). The cross section does not
max forcecanassumedifferentexpressions,dependingonthephysical
varybymuchwhenthedistributionistruncatedtosentr (cid:38) λ,be- regime.Ifthesizeofthedustparticles(cid:46)λ ,themean-freepath
mfp
causethelargestcontributiontothecrosssectioninthegeometric
ofgasmoleculeswithintheflow,thedragforcereducestotheso-
limitcomesfromthesmallestgrainswhenq > 3.Wedonotre-
calledEpstein(1924)drag.Ifthesizeoftheparticleislargerthan
computethetemperatureT(n)withthePDRcodeforthenewtrun-
the mean-free path, drag is a classical fluid Stokes (1851) drag.
cateddistributions,sinceitdependsweaklyonthemaximumgrain Withinthewindsaddressedinthispaper,λ > 105cm(using
mfp
size.AsdiscussedinSection3,thisisduetothefactthatheatingis λ ≈ (σ n)−1,whereσ ∼ 10−16cm2 isthegeometrical
mfp gas gas
mostlygeneratedbythephotoelectriceffectonPAHsandonsmall cross-section of the gas moleculesand n < 1011cm−3). For the
grains,sincethetotalsurfaceareaperunitmassisdominatedby
dustparticlesconsideredhere,s(cid:28)λ ,andwecanthereforeuse
mfp
the small grains of the distribution when q > 3. Cooling is not
theEpsteinlimittoestimatethedragterm(e.g.Armitage2010):
significantlyaffectedbytheabsenceofthelargestgrains,sinceOI
andCIIemissionlinesdominatetheradiativecooling. F = 4πρs2v (v −v), (29)
D 3 th j
whereρ=µmnisthegasmassdensity,and
5.2 Fluidequationsfordustparticles
(cid:115)
Inordertoobtainthemaximumgrainsizeentrainedinagaseous v = 8kBT (30)
solution, we need to solve the hydrodynamical equations of dust th πµmH
(cid:13)c 2015RAS,MNRAS000,1–18