星形成ゼミ
2018/06/08
SFN #304 #49-#58
長谷川哲夫（国立天文台）

1

Scottish Universities Physics Alliance (SUPA), School of Physics and Astronomy, University of St. Andrews, North Haugh, St.
Andrews, Fife KY16 9SS, UK

304 Mul$-epoch monitoring of the AA Tau like star V 354 Mon#49 Indica$ons for a low gas-to-dust ra$o in the inner disk warp

received; accepted

P. C. Schneider1; 2, C. F Manara2; 3, S. Facchini4, H. M. Günther5, G. J. Herczeg6, D. Fedele7, and P. S. Teixeira8; 9; 10
A&Ap in press
ABSTRACT

A&A proofs: manuscript no. final
Disk warps around classical T Tauri stars (CTTSs) can periodically obscure the central star for some viewing geometries. For these soTau-like variables,
obscuringis located in the inner disk and absorption spectroscopy allows one to characterize
arecalled
not AA
necessarilyaligned the(Ansdell
et materialal.
its dust
andabsorption
gas content.by
Since
thediskobserved
emission from CTTSs consists of several components (photospheric, accretion, jet, and
2016)
and
inner
matedisk emission), which can all vary with time, it is generally challenging to disentangling disk features from emission variability. Multiin epoch,systemswith rather inclined outer disks,
flux-calibrated, broadband spectra provide us with the necessary information to cleanly separate absorption from emission
ypicalAA TauWesystem,
ALMA
obser-epochs of VLT/X-Shooter spectra of the CTTS V 354 Mon (CSI Mon-660) located invariability.
appliedwhere
this method
to threeouterreveala diskof about
NGCdisk
2264
and findthat:inclination(a) the accretion
emission remains virtually unchanged between the three epochs; (b) the broadband flux
re evolution
decidedlyis not
to edge-on(Loomisbest close
described
by disk material
obscuring part of the star, and (c) the Na and K gas absorption lines show only a minor
increase in equivalent width during phases of high dust extinction. The limits on the absorbing gas column densities indicate a low
ratio in no
the strong
inner disk,
lessthan a tenth of the ISM value. We speculate that the evolutionary state of V 354 Mon, rather
he gas-to-dust
AA Tau system,colorchanges
withfound.
a low accretionis responsible
aveoldbeen
In other rate,words,
the extinc-for the dust excess through an evolution toward a dust dominated disk or through thefragmentationof larger
that drifted
inward from larger radii in a still gas dominated disk.
whichsuggests thatthe bodies
extinction
is caused

e linesight Stars:(Bouvier
et al. V354
2003).
HowKey ofwords.
individual:
Mon,
Stars: low-mass, stars: pre-main sequence, stars: variables: T Tauri, Herbig Ae/Be, X-rays:
bility can also change the observed stellar.k-locking (e.g., Matt et al. 2010), the stellar-on period of the inner disk warp are likely
変光周期 5.25日,
innerdisk structure
CTTSs indicating
ourtechniques. Tedsothat it is challenging to disentangle ac- Fig. 1. Sketch of theIntroduction
rentlyunaccessible
byaround(interferometric)
imaging
threecomponent
model.
The
unocculted
part
of
the
stellardisk
isdent properties using photometry
only.
periods of the dips are typically five to ten days,
which translaiA surrounded by copious amounts of notedwith A. The upper disk layer causes reddening of part the stellar
s
ung
stars
are
circumstellarhwarps are not stable and changes in the inner
radiidisk).
of about
au assuming
Keplerianrotation. Hence,(part
Bacof the tostellar
Stellar0.1
emissionfrom theregion denoted
B fraction of this material survives disk
and
dust.
A
the
main
se preferred explanation for major dimming
warp
material
close to byor theeven
cospatialwith the inner dwithC isa subject
to opaque(gray)isextinction
lower
disk layers.
C can still surround the young star
tion
phase
and
after
few
Figure1 sketches
the envisaged
geometry,
damples of theser dimming events are AA Tau
We
also mightshow the edge.
inner disk
radiusr, the viewing
inclinationi, and thewhich natura
lion
years
in
the
form
of
a
circumstellar
disk,
which
-; Schneider et al. 2015a), RW Aur (Petrov
results
fromwiththe respect
interaction
the stellar magnetic fiheight h of the innerdisk warp
to thebetween
disk midplane.
~0.1 AU
rlve into a debris
diskdominatedbyreprocessed
dust
grains.
and the inner disk edge (Foucart & Lai 2011).
er
et al. 2015b; Bozhinova et al. 2016; Faceyr wellstructureofprotoplanetarydisks
iscurrently
subject
toina few other objects (Hamilton et al.
1. Sketch of the innerdisk
structure around
CTTSs indicating our often focusingseobservational
andtheoreticalinvestigations,- Fig.
A related, perhaps overlapping, class of objects are the(a)accreQon
signaturesl. 2015,2017).
Combiningthefindings
threecomponent
model.
The
unocculted
partof
the stellar fromdiskusingis de- emission line ﬂuxes
成分分離
intrinsic
emission(Nomuraetofal.
Ansdellettheal.broadband
called dippers(Stau↵er
et al.
2015). (SED)
Theirmaterial,
light curves bear maulationinof
theobserved
stellarflux
due tocircumstellar
noted with disk
A. The upperdisk layer
causesreddeningpart 2016;
the stellar(b)dust
exQncQon
throughchanges
spectralenergydistribuQon
r use
information
from
X-ray
tofromNIR
wavedisk
(part
B
of
the
stellardisk).
Stellaremissiontheregiondenoted16;Cox
et al. 2017;
Loomisetcaused
al.
2017;
Long
ettoal.obtain2017).further
similarities
to thethose
AA Tau-likebut the dipsare lg withinsight
intocompositionof theobjects,
disk
C is subjectevents
to opaque(gray)
extinction
by thethe lower
layers. gas
2(c)arechangesinabsorbingcolumndensityfrom inspecQon
ofofstrong
absorpQonlinesmatehe dimming
likelybydiskperu AnWe alternative
also show the inner
radius r,disk
the viewing
inclinationstudies
i, and the is absorption
todisktheseemission
pronounced and shallower than for those of “classical” AA Tera scaleheightof about
stel- rial.
v diskheight hwithof the innerdisk warp
with respect
to the diskone
midplane.

manuscript
finalstars,no.Infrared:stars, Protoplanetary disks観測： VLT X-Shooter（λ300 – 2500 nm分光測光）、連続3夜
結果：1) 輝線（Hα、HeI、Li）弱い→accreQonの変化ではなくoccultaQonの変化
2) SEDは大きく変化→warped inner diskのdustによる減光
3) 吸収線（NaとK）→大部分photosphericだがinner diskによる変化もある
Gas/dust比はinterstellarに比べて1桁以上低い
仮説1: inner diskはdebris diskに進化中
P. C. Schneider et al.: Multi-epoch monitoring of the AA Tau like star V 354 Mon
仮説2: outer diskから落ちてきた粒径の大きな個体の破砕

Fig. 3. Top: Observed spectral energy distributions. Bottom two panels: Absolute di↵erence between model and observed spectra. Models for thedimmer spectra (18th and 19th Jan.) assume that the di↵erence between the brightest spectrum (20th Jan.) and the dimmer spectra is caused (1)only by additional extinction allowing E(B V) and RV to vary (blue), and (2) our AA Tau model in red (see text for a detailed description). To
di↵erentiate between the two 18th Jan. spectra, the shorter one is displayed as a dashed line. Spectral ranges strongly a↵ected by telluric absorption
are removed from the fit and shown in lighter colors. The parameters corresponding to the blue and red curves in the bottom two panels are given
in Table 2.

3

Patrick D. Sheehan

and Josh A. Eisner

Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

304 MulQple
Gaps
inJulythe11;Diskof2018
theFebruary
Class I7;Protostar
GYFebruary
91 7; published 2018 April 9
Received
2017
revisedaccepted 2018
D. Sheehan and Josh A. Eisner#50 Patrick
ApJ 857, 18
Abstract

We present the highest spatial resolution Atacama Large Millimeter/submillimeter Array (ALMA) observations todate of the Class I protostar GY 91 in the ρ Ophiuchus L1688 molecular cloud complex. Our 870 μm and 3 mmdust continuum maps show that the GY 91 disk has a radius of ∼80 au, and an inclination of ∼40°, but mostinterestingly that the disk has three dark lanes located at 10, 40, and 70 au. We model these features assuming they
The Astrophysical Journal, 857:18 (7pp), 2018 April 10& Eisneraregaps in the disk surface density proﬁle and ﬁnd that their widths are 7, 30, Sheehan
and 10au. These gaps bear a striking
resemblance to the gaps seen in the HL Tau disk, suggesting that there may be Saturn-mass planets hiding in the
disk. To constrain the relative ages of GY 91 and HL Tau, we also model the disk and envelope of HL Tau and ﬁndthat they are of similar ages, although GY 91may be younger. Although snow lines and magnetic dead zones canalso produce dark lanes, if planets are indeed carving these gaps then Saturn-mass planets must form within theﬁrst ∼0.5 Myr of the lifetime of protoplanetary disks.
Key words: planets and satellites: formation – protoplanetary disks – stars: formation – stars: individual (GY 91)(n5–12=−0.25). However, their measured value is athe border between disks with foreground extinction andnets form in protoplanetary disks. When observed at highwith envelopes (of n5–12=−0.2), and the extinction-cor
ion with Atacama Large Millimeter/submillimeter Array
spectral index (a ¢IR =0.31) and bolometric temp
A),Figurea number
these
disks
show
interesting
patterns
1. Our ALMAof
345 GHz(left) and
100 GHz(right) mapsof the GY 91 protoplanetary
disk. inTwo dark lanes are readily apparent in the 345 GHz map, while a¢ =470(TGHzK) still qualify the source as a Class I pr
third dark lane is also apparent in inner regions of the disk at 100 GHz because of the higher resolution of our 100
maps.
bol
millimeter emission proﬁles, in some cases including(Dunham et al. 2015). van Kempen et al. (2009) alsoof bright and dark rings (ALMA PartnershipTable
et 1al. 2015;
HCO+ emission toward GY 91 that was bright enoughLog of ALMA Observations
et al. 2015; Andrews et al. 2016; Isella et al. 2016; Fedele
cutoff for an embedded source, but that em
Observation Date
ALMA Band
Baselines
Integration Time above theCalibrators
2017a,
2017b; Loomis
et al. 2017; Dipierro
et al. Total2018).
(UT)(m)(s)(Flux,
Bandpass,
Gain)
seems to be associatedwith a patch of cloud that peaof these
disks
may
even
be
quite
young
(1
Myr;
ALMA
2015 Oct. 31
3
84–16, 200
169
1517–2422, 1625–2527
away fromthe source.
2015 Nov. 26
3
68–14, 300
169
1517–2422, 1625–2527
rship
et
al.
2015;
Dong
et
al.
2015;
Fedeleetal.
2017b;
2016 Apr. 17
3
15–60058
1733–1304,
1427–4206,new1625–2527
HereJ1517–2422,
we present
ALMA data, which, when com
May 19730
J1625–2527
ro 2016et al.
2018). There are
a number15–640
of explanations forwith theJ1517–2422,
observed
SED, show that GY 91 is indeed a C
2016 Sep. 11
7
15–314060
J1625–2527
eatures, including chemical processes that alter dust
protostar with a circumstellar disk embedded in4 an env
es wereanddevoted
sticking/fracturing
processes
near snow
et al. lines
2002; Natta et al.
2006;
Stankeandet al.870
2006;μm
Alves
deto spectral line observations,although nothing
Our
3 mmimagesalso reveal the presence oOliveira & Casali 2008; Jørgensen et al. 2008; Padgett
was detected. We list details of the observations in Table 1.
1. Introduction

The Astrophysical Journal, 857:18 (7pp), 2018 April 10

Sheehan & Eisner
The Astrophysical Journal, 857:18 (7pp), 2018 April 10

Figure 1. Our ALMA 345 GHz (left) and 100 GHz (right) maps of the GY 91 protoplanetary disk. Two dark lanes are readily apparent in the 345 GHz map, while a Figure 2. One-dimensional, azimuthally averaged, de-projected radial brightThe dark
Astrophysical
April
Sheehan & Eisner ness proﬁle of the GY 91 disk at 345 and 100 GHz, with the locations of thethirdlane is also Journal,
apparent in857:18
inner (7pp),
regions2018of thedisk10at 100 GHz because of the higher resolution of our 100 GHz maps.

dark lanes marked by vertical dashed lines. These gaps are readily seen in thebrightness proﬁle. We also show the azimuthally averaged brightness proﬁle ofour gapped disk+envelope model (see Figure 4, Table 2) at each wavelength.

Table 1
Log of ALMA Observations
Observation Date
(UT)2015 Oct. 31
2015 Nov. 26
2016 Apr. 17
2016 May 19
2016 Sep. 11

ALMA Band

Baselines
(m)

3
3
3
7
7

84–16, 200
68–14, 300
15–600
15–640
15–3140

bins. The outer two dark lanes show up clearly in the 870 μm
radial proﬁle, while the inner lane shows up clearly in the 3 mm
proﬁle. Moreover, there appears to be a break in the 870 μm
proﬁle at the location of the inner dark lane, and there appears
169
1517–2422, 1625–2527
to be a dip in the 3 mm brightness proﬁle that is consistent with
169
1517–2422, 1625–2527
the location of the middle dark lane, despite the noisiness of the58
1733–1304, 1427–4206, 1625–2527 3 mm image that prevents it from being detected by eye.
30
J1517–2422, J1625–2527
In order to study these features in greater detail, we ﬁt a
60
J1517–2422, J1625–2527
model to the data to determine disk properties such as radius,
position angle, and inclination, as well as the locations, widths,
and depths of the gaps. We use Monte Carlo radiative transfer
et al. 2002; Natta et al. 2006; Stanke et al. 2006; Alves de codes to produce synthetic observations of model protostars
Oliveira & Casali 2008; Jørgensen et al. 2008; Padgett that can be ﬁt to our combined millimeter visibility andet al. 2008; Wilking et al. 2008; Evans et al. 2009; Gutermuth broadband SED data set of GY 91. This modeling procedure is
et al. 2009; Barsony et al. 2012). When modeling the SED, as described in further detail in Sheehan & Eisner (2014) and
we discuss below, we assume a constant 10% uncertainty on Sheehan & Eisner (2017), but we give a brief overview here.
Our model includes a ﬂared protoplanetary disk with aany photometry from the literature to account for any ﬂux physically motivated surface density proﬁle (e.g., Lynden-Bellcalibration uncertainties between the measurements.
& Pringle 1974) surrounded by a rotating collapsing envelope
In addition to the broadband photometry, we also download (e.g., Ulrich 1976),
Total Integration Time(s)

Calibrators(Flux, Bandpass, Gain)were devoted to spectral line observations, although nothing
was detected. We list details of the observations in Table 1.
The data were reduced in the standard way with the CASA
pipeline and the calibrators listed in Table 1. After calibrating, we
imaged the data by Fourier transforming the visibilities with the
CLEAN routine. After our initial imaging, we found that we could
improve the sensitivity of the 345 GHz image by self-calibrating.
We ran four iterations of phase-only self-calibration on thethe Spitzer IRS spectrum of GY 91 from the CASSIS database
compact conﬁguration track and a single iteration of phase-only(Lebouteiller et al. 2011, 2015). Rather than consider the entire
self-calibration on the extended conﬁguration track. This
SED, which can be computationally prohibitive for the radiative
improved
the rms in an image produced with natural weighting
Figure 3. Best-ﬁt simple geometrical model for GY 91 compared with−1
the data. The model assumes the disk is ﬂat, with a surface density described by Equation (1)transfer
calculations
below,(i.e.,
robustparameterof au,
2) γ=0.3,
from 0.36
to 0.27
beam The. Wemodel uses
T = 46we(R 1sampleau)-0.4 . Wetheand Madisk
=0.36
Me, rc=71
i=39°,and mJy
p.a.=−19°.
a power-law
temperature described
distribution withuse IRSthe2 −1at)b25
ranging
5 to 35
μm. three
We gapsalso withassumewereunableto improve
100described
GHz imageby self-calibration.
cmpoints
g withβ=1.8.from
Our model
includesthepower-law
millimeter
opacity thefunctionin Beckwith
et al. (1990), k (n ) = 0.1 (spectrum
n 1000 GHz
au, usingwgap,1=5.9
au, weightingδgap,1≈0, with
Rgap,2=40.3
au,
wgap,2uncertainty
=27.5 au, δgap,2Rgap,3=68.9
au, andfollowingparameters:
Rgap,1=10.4
a 10%
on =0.15,
these ﬂuxes,
like weau,dowgap,3for =10.7
the broadband
Our ﬁnal
images were
produced
Briggsa
Our modelingindicates
thatprovides
the ﬁrst andthird gaps
are deep;
however, as thedata is noisy and does not have a high enough resolution to well resolve theδgap,3≈0.
photometry.
robustparameter
of
0.5,
which
a
good
balancebetweengaps, the actual depths are quite uncertain. As the disk is quite massive, it is somewhat optically thick at 870 μm. We show the one-dimensional, azimuthally averaged,
sensitivity
and resolution,
theimages
visibilities
for both thedatavisibility amplitudes
on the left,totheweight
raw modelon the center-left,
model images convolved with the beam in the center-right column, and the residuals on
sets.
The The3 mm
image
hasthea beamof sizeby 0 05andwitha P. panels. The peak residuals are 1.7σ at 345 GHz and 3.5σ at 100 GHz.
the right.
colorbarshowscolor scalefor 0
the 06center-right
right-most−1
3. Results
A. of 81°. 9 and an rms of 36 μJy beam . The 870 μm image has⎡ ⎛ R ⎞2 - g ⎤⎛ R ⎞-g
S = S0 ⎜ ⎟ exp ⎢ - ⎜ ⎟ ⎥ ,⎢⎣ ⎝ rc ⎠ ⎥⎦⎝ rc ⎠(1 )⎛ 1 ⎡ z ⎤2 ⎞S
exp ⎜ - ⎢ ⎥ ⎟ ,⎝ 2 ⎣h⎦ ⎠
2p h

(2)

rdisk =

⎛ R ⎞b⎟ ,h = h0 ⎜⎝ 1 au ⎠
1

5

(3 )

The Astrophysical Journal, 857:18 (7pp), 2018 April 10
The Astrophysical Journal, 857:18 (7pp), 2018 April 10

Sheehan & Eisner
our gapped disk+envelope model (see Figure 4, Table 2) at eachwavelength.
Sheehan& Eisnerbins. The outer two dark lanes show up clearly in the 870 μm
radial proﬁle, while the inner lane shows up clearly in the 3 mm
proﬁle. Moreover, there appears to be a break in the 870 μm
proﬁle at the location of the inner dark lane, and there appears
to be a dip in the 3 mm brightness proﬁle that is consistent withthe location of the middle dark lane, despite the noisiness of the
3 mm image that prevents it from being detected by eye.
In order to study these features in greater detail, we ﬁt a
model to the data to determine disk properties such as radius,
position angle, and inclination, as well as the locations, widths,
and depths of the gaps. We use Monte Carlo radiative transfer
codes to produce synthetic observations of model protostars
that can be ﬁt to our combined millimeter visibility andbroadband SED data set of GY 91. This modeling procedure is
Figure 4. Gapped disk+envelope model for GY 91 compared with the data. We show the one-dimensional, azimuthally averaged, 870 μm visibility amplitudes on the
Figuredisk+envelopemodel forcomparedthe one-dimensional,
averaged,
870 μmthedescribed
in furtherazimuthally
detail in
Sheehan&visibility
Eisneramplitudes(2014) onand
left, the4.3Gapped
mm visibilities
in the center,
andGY
the91
SED
on the withright.the data. We show
left, the 3 mm visibilities in the center, and the SED on the right.
Sheehan & Eisner (2017), but we give a brief overview here.
Our model includes a ﬂared protoplanetary disk with a
Johansen 2013; Banzatti et al. 2015). The change in the optical
Table 2
physicallymotivatedsurfacedensity
proﬁle(e.g.,
Johansen2013;
Banzatti
et al.
2015).
The
changeincould
the optical
Table 2 Model Parameters
Gapped Disk+Envelope
propertiesof dustgrainsacrossthesnowlineLynden-Bell
cause
Gapped Disk+Envelope Model Parameters& Pringle
1974)
surroundedbyarotating
collapsingenvelope
properties
ofdust
grainsacross
the
snow
line
could
cause
features like those seen in HL Tau or GY 91 (e.g.,
Zhang
Parameters
Valuesfeatureslikethose seen in HL Tau or GY 91 (e.g., Zhang(e.g.,
1976),
Parameters
Valueset Ulrichal. 2015).
Lstar [L :]0.11±0.00
et al.
2015).
Alternatively, the “sintering”
dust
just below the
2 -grainsg⎤LM
0.11±0.00⎡ ofof⎛dust:] ]star [L[M⎞-gproduces
0.17±0.00
Alternatively,the⎛ R“sintering”
just
the
disk
e
R ⎞ grains
sublimation
temperature
brittlegrainsthatbelow
fragment⎥,[Me]
0.17±0.00
M(1 )S temperature= S0 ⎜ ⎟ produces
exp ⎢ - ⎜brittle⎟ grains
0.1±0.0
Rdisk
in [au]sublimationthatfragment
more readily and therefore
smaller⎝ rc ⎠ growrc ⎠ ⎦⎥sizes just outside of
0.1±0.0
RRin [au]⎣⎢ toto⎝smaller
75±1disk [au]morereadily
andBecause
therefore
outsideofthe snowline.
thegrow
sinteredgrainssizes
havejust
smallersizes,[au]
75±1
Rhdisk[au]0.29±0.01
0
the
snow
line.
Becausethesinteredgrainshave
smallersizes,
0.29±0.01
hγ0 [au]they undergo slower Sradial drift,
pile-ups near snow
⎛ causing
0.01±0.04
1 ⎡ z ⎤2 ⎞pile-ups
theyslower radialexp
drift,
neartosnow
γβ0.01±0.04(2)
lines.undergo
Thisrprocess
producefeaturesthose⎟ , similar⎜ -causing
0.30±0.01⎢⎥
disk = could also⎣⎦lines.
This
process
could
also
producefeaturessimilar
2
h
⎠
⎝
p
2
h
βM [M ]
0.30±0.01
seen in HL Tau or GY 91 (Okuzumi et al. 2016). to those
0.020±0.008
enve
seenWein compare
HL Tau theor GY
91 (Okuzumi
et al. 2016).
M[Me]
0.020±0.0081269±328
Renv
midplanedisk btemperatureinferred for GY
env [au]1269±328
Rfcav
We
comparethe
midplanedisktemperatureinferred volatiles
for GY⎛⎞env [au]
R
0.56±0.04
91 with the temperatures
of
snow
lines
of
common⎜⎟ ,(3 )h = h 0of
0.56±0.04
fcav
91
with
the
temperatures
snow
lines
of
commonvolatilesξ
1.31±0.24(e.g., Zhang et al. 2015).⎝ 1The
au ⎠outermost gap does roughlyξi [◦]
1.31±0.24(e.g.,
et al. 2015).
The ofoutermost
gap does
41±1
matchZhang
the freeze-out
regionN2. However,no roughly
obvious[◦] [◦]
ip.a.
41±1
1 . However, no-1obvious
match
the freeze-out
regionofN
111±1
2
counterparts
are
seen
for
the
inner
two
gaps.
2
◦
⎛
⎛m
111±1p.a.
m ⎞ twoM˙ are seen
R ⎞82597±6318
amax[ ][μm]counterpartsinner3)- 12for
ﬂows(GM*rproduced
renvZonal=+ by+ 2m20 c ⎟ instability. (4 )⎜1 the⎟magneto-rotational⎜ gaps.
[μm]
82597±6318
aRmax
[au]
10.0±0.2
gap,1
Zonalturbulencebyinstabilitym
4pﬂows producedr ⎠been⎝⎠magneto-rotational⎝ m 0 have also
driven(Johansen
et0 al.
2009)
shown
Rgap,1 [au]10.0±0.2
wgap,1 [au]7.0±0.4
driventurbulence(Johansenpressure
et al. 2009)
havethat
alsocanbeen
shown
to produce
axisymmetric
bumps
trap
large
wgap,1 [au]7.0±0.4
dgap,10.01±0.02 In Equations
todustproduceaxisymmetricthat cantrap
large(1)–(3),
R produce
andpressurez aregap-likeinbumpscylindricalcoordinates,
grainsandmay
featuresinmillimeter
dgap,10.01±0.02
dust
and (4),may
produce
gap-likefeaturesinin millimeter
Rgap,2 [au]40.9±0.4
while
ingrains
Equationμ=cos
and r and
θ arespherical
images(Pinilla
et al.
2012; θDittrich
et
al.
2013;
Simon &
Rgap,2 [au]40.9±0.4
images(Pinilla
et
al.
2012;
Dittrich
et
al.
2013;
Simon
&
wgap,2 [au]30.0±6.7
Armitage2014).
In
this
case,
the
pressurebumps
arecreatedbycoordinates.
wgap,2 [au]30.0±6.7
Armitage2014).
In thisincase,the pressure bumps
are created thatbydgap,20.16±0.07
large-scalevariations
the magneticallydriven turbulence
dgap,20.16±0.07
large-scalevariationsininthethemassmagneticallythat 3
Rgap,3 [au]69.4±0.6
produce variations
accretion driven
rate thatturbulencein turn causes
Rgap,3 [au]69.4±0.6
producevariations
in
the
mass
accretion
rate
that
in
turn
causeswgap,3 [au]10.0±0.8
material to pile up. This effect can also be seen at the outer
wgap,3 [au]10.0±0.8
material
pile up. deadThis effectalso there
be seenat the outeredge of tomagnetic
zones,canwhereis strongradialdgap,30.01±0.02
edgeofmagnetic
deadzones,
where
there
is
strongradialdgap,30.01±0.02
variation
in
the
mass
accretion
rate.
These
ﬂows
can produce
p3.50
variationin the massaccretion
These
producep
3.50
gap-like featuresin disks(Pinillarate.
et al.
2012;ﬂows
Flockcanet al.
2015).
gap-like features in disks (Pinilla et al. 2012; Flock et al. 2015).

Disk+envelopeのモデル

as initial guesses to generate a detailethat simultaneously matches both vi
broadband SED.
The purpose of the radiative transf
is to determine whether an envelope
explain the observations. We use thtransfer codes RADMC-3D (Dullem(Robitaille 2011) to calculate the tedensity structure described above,
synthetic millimeter visibilities andthese synthetic observations simultan
data sets (870 μm visibilities, 3 mm v
SED), again using the emcee code.
The MCMC walkers are allowed tospace for an extended burn-in perio
have converged when the walkers hawith the best-ﬁt value for each paramover a large number of steps. We t
value for each parameter as the media
after discarding the burn-in steps. W
on these parameter values as the rcontaining 68% of the walkers.
We compare the gapped disk+
observed visibilities and SED in Fig
parameters in Table 2. The imag
+envelope model look almost iden
Figure 3, although the residuals betw
higher. This is not surprising, as weand SED simultaneously here.
This model can simultaneously
visibilities, 3 mm visibilities, and bro
The GY 91 disk appears to be embe
Menv=0.12 Mdisk. Gaps are found aand ∼69 au, with widths of ∼7 au, ∼gap depths for the inner and outer gap
because they are not well-resolved
middle gap appears to be wide
although with higher resolution it is
up into multiple gaps.

6

ESPOL Polytechnic University, Escuela Superior Politécnica del Litoral (ESPOL), Physics Department, FCNM, Campus Gustavo Galindo km. 30.5
Vía Perimetral, P.O. Box 09-01-5863, Guayaquil, Ecuador
Received2018 March 1; accepted2018 March 8; published
AprilBinary
11
304 From Large-scale to Protostellar
Disk FragmentaQon
into2018
Close
Stars
Di G. Sigaloi1, Fidel Cruz1 , Ruslan Gabbasov2, Jaime Klapp3,4 , and José Ramírez-Velasquez5,6#51 Leonardo
ApJ 857, 40
Abstract

Recent observations of young stellar systems with the Atacama Large Millimeter/submillimeter Array (ALMA)
and the Karl G. Jansky Very Large Array are helping to cement the idea that close companion stars form via
fragmentation of a gravitationally unstable disk around a protostar early in the star formation process. As the disk
grows in mass, it eventually becomes gravitationally unstable and fragments, forming one or more new protostars
in orbit with the ﬁrst at mean separations of 100 au or even less. Here, we report direct numerical calculations down
to scales as small as ∼0.1 au, using a consistent Smoothed Particle Hydrodynamics code, that show the large-scale
fragmentation of a cloud core into two protostars accompanied by small-scale fragmentation of their circumstellar
disks. Our results demonstrate the two dominant mechanisms of star formation, where the disk forming around a
protostar (which in turn results from the large-scale fragmentation of the cloud core) undergoes eccentric (m=1)
fragmentation to produce a close binary. We generate two-dimensional emission maps and simulated ALMA
1.3 mm continuum images of the structure and fragmentation of the disks that can help explain the dynamical
processes occurring within collapsing cloud cores.
Key words: circumstellar matter – methods: numerical – planets and satellites: formation – protoplanetary disks –stars: formation – stars: pre-main sequence
1. Introduction

lar systems, consisting of two or more stars, are formedo different mechanisms—large-scale fragmentation ofcores of gas and dust during their early isothermal
se phase (Larson 2001; Fisher 2004), leading to widely
ted companion protostars (∼1000 au or larger), andscale fragmentation of their circumstellar disks early inollapse due to gravitational instabilities (Rodríguez
2005; Connelley et al. 2008; Pech et al. 2010; Murillo
2013; Tobin et al. 2013), leading to bound systems
ted by tens of au. This scenario of star formation seems

lack physical reality because protostellar disks form andin mass within a collapsing environment as a result of theof material from the outer cloud envelope (Eisner 2012;
et al. 2012). On the other hand, if more protostars form w
the same environment, the gravitational interaction athem may limit the disk size and affect the way anmomentum is transferred. However, the disk fragmen
model has succeeded in reproducing the constraints impos
the observed statistical properties of low-mass objects, su
low-mass binary systems, brown dwarfs, and even planmass objects, which in turn have not yet been7 explain
means of other formation mechanisms (Stamatellos & Whi

The Astrophysical Journal, 857:40 (11pp), 2018 April 10

Sigalotti et al.

Figure 1. Equatorial density evolution in a 266 au region around the cloud center for model U2 after large-scale fragmentation into a binary protostar: (a) 44.14 kyr,(b) 44.39 kyr, (c) 44.43 kyr, and (d) 44.68 kyr. The logarithm of the density ρ is plotted. For (a)–(d), the maximum density is: log10ρmax=−10.86, −10.63, −10.62,
and −10.39 (ρmax in g cm−3), respectively. The insets show a magniﬁcation of the binary cores and their whirling disks. The disks develop a two-armed spiral patternassociated with a linear growth of a gravitational instability. In (b) and (c), the gravitational instability in disk B enters a nonlinear growth phase and so one of its spiralarms is seen to fragment into a second protostar.

8

Figure 2. Time variation of Toomre’s Q parameter along the radial structure of disks A (left) and B (right). The color bar in the middle of both plots indicates thevalues of Q. Disk A shows values of Q>0.6 and hence it does not fragment. In contrast, Q<0.6 is achieved in one of the spiral arms of disk B as highlighted by the
undulating yellow region between 10 and 20 au from the primary. This agrees with a recent revision of the Toomre’s criterion, which indicates that Q<6 in the spiralarms is a necessary condition for fragmentation of the disk (Takahashi et al. 2016). As a result, the disk undergoes eccentric (m=1) fragmentation into a second
protostar, forming a close binary.

Figure 3. (a) Time evolution of the masses of the primary protostars (thin solid blue and dashed green lines) and disks (thick solid blue and dashed green lines)
produced in model U2 for systems A and B. (b) Disk mass over central protostar mass ratio, Md/Må, for disks A and B. The thick solid red line in (a) depicts the mass
involved in the gravitationally unstable region within disk B as a function of time. The horizontal dotted line in (b) marks the analytically derived limit Md M = 0.3
above which eccentric (m=1) disk fragmentation is expected to occur (Adams et al. 1989). Disk B grows to Md/Må≈0.7, while disk A exhibits values of Md/Må
that are barely above the limit for eccentric fragmentation.

9

The Astrophysical Journal, 857:40 (11pp), 2018 April 10

Sigalotti et al.

The Astrophysical Journal, 857:40 (11pp), 2018 April 10

Sigalotti et al.

Figure 4. (a) Projected particle positions on the f=0 plane of a cylindrical pillbox of radius 50 au and height 40 au enclosing disk A. The vertical red lines binds
three bins of radial width 1 au at selected radii from the disk center. (b)–(d) Gaussian-like density distributions corresponding to the bins in (a) as indicated by thecolored arrows. (e) The same as (a) but for disk B. (f)–(h) The same as (b)–(d) for disk B. The maximum and width of the distributions deﬁne the disk midplane(R) and scale Journal,height H(R)usedin Equation(3),
densityρ0AstrophysicalThe
857:40(11pp),
2018 April
10 respectively.
Sigalotti et al.

we calculate the temperature for all particles by combining thebarotropic pressure–density relation with the ideal gas law to givecollapse was seen to remain strictly isothermal up toρ∼10−16 g cm−3. At higher densities, the collapse was nearisothermal, with the temperature rising very slowly from 10 K
g1
⎡
⎤
⎛r⎞
Particle density distribution as a function of radius for (a) disk A and (b) disk B. The horizontal solid line marks the transition from isothermal tohigher5.thanto ∼11.2 K when ρ∼10−14 g cm−3. At densities Figure
T = T0 ⎢1 + g ⎜ ⎟ ⎥ ,(4 )nonisothermal collapse. Also shown is the particle temperature distribution as a function of radius for (c) disk A and (d) disk B. The time shown in all plots is⎢⎣
this, the temperature was seen to rise sharply, witht=44.39
the density⎝ rc ⎠ ⎥⎦kyr. The density and temperature peaks centered at R=0 represent the primary protostars, while those centered at R≈18.2 au correspond to the second
fragment
producedin disk B.
and the collapse becoming nonisothermal. Therefore,a valueofwhere T0=10 K, rc = 5.0 ´ 10-12 g cm−3, and γ=5/3. At
the critical density rc = 1.0 ´ 10-14 g cm−3 would be moreto each cell center by simply averaging the density of all of
thosea valueat R≈18.2 au correspond to the secondary fragment
temperatures below 100 K, a value of γ=5/3 is appropriate
representative of the near isothermal phase. However,those particles within a cell. This procedure permits drawing
in disk B. The upper portions of the density peaks arebecause the rotational and vibrational degrees of freedom ofisothermalof rc = 5.0 ´ 10-12 g cm−3, while prolonging theproduced
density proﬁles from cellwise densities along the z-direction
well above the line ρ /ρc=1, implying that the protostars aremolecular hydrogen are frozen out, and so only translational
phase of collapse to higher densities, makes convergenceover the entire volume of the box, which are then ﬁtted to a
in a nonisothermal phase of collapse. The small peaks present
degrees of freedom need be considered (Boss et al. 2000).
more difﬁcult to attain and resolution requirements signiﬁcantly
Gaussian distribution to determine the disk midplane density,about
10 au from the central protostar in Figures 5(a) and (c)However, precise knowledge of the dependence of temperature
more demanding. Larger and more massive disksatthanthoseρ0(x, y), and scale height, H(x, y) (see Section 3.1). The disk
clearly indicatethat by this time, disk A is in the process of
on density at the transition from isothermal to nonisothermal
calculated here can therefore be obtained by anticipatingthesurface density Σ (in g cm−2) is then calculated from thefragmenting
collapse will require solving the radiative transfer problem
near isothermal phase using lower transition densities
in into a new secondary, as can be better seen fromestimates of H(x, y) using the expression Σ=2ρ H, where ρ isemission maps of Figure 6(a).
coupled to a fully self-consistent energy equation. Thisbarotropic collapse calculations. Figure 5 shows thethe
projection
the gas volume density (Dong et al. 2016). The disk midplane
The images displayed in Figure 6 correspond to disk surface
transition was studied for the collapse of an initially centrally
temperature is ﬁnally calculated using Equation (4) with ρon the disk midplane of the disk particle densities (normalized
density and midplane temperature maps at selected times
condensed core obeying a spherically symmetric Gaussian
replaced by ρ0(x, y) as determined from the maximum of theto the critical density rc = 5 ´ 10-12 g cm−3) and temperaturesduring the evolution of disks A and B. These images are
Gaussian z-density proﬁles.
density proﬁle, using nonisothermal thermodynamics and(calculated using Equation (4)) as functions of radius
for disks
generatedfrom the numerical data by ﬁrst drawing a cubic box
The radiative ﬂux emerging from the disk surface is ﬁnallysolving the mean intensity equation in the Eddington approxA and B at ≈44.39 kyr. The prominent peaks in density
andof length 56 au that encloses the entire disk, wherewith sides
calculated using the following relation (Dong et al. 2016):
imation with detailed equations of state (Boss et al. 2000). The
temperature at R=0 represent the central protostars,
while
the geometricalcenter (x=y=z=0) and the z=0 plane of
the box are made to coincide with the center-of-mass of the
4tP sT 4,(5 )
F (S1 2 ) =primary core and the approximate equatorial plane of the disk,
3
6
1
2
+
tP + 2 tR tP
respectively. The box is divided into a number of parallelz=const. planes above and below the z=0 plane, with
Figure 8. (a) Time evolution of the masses of the primary Må (solid blue line) and secondary MGI (solid magenta line) protostars in disk B. (b) Orbital evolution of thewhere σ is the Stefan–Boltzmann constant, T is the disk
separations of 0.35 au between each other. Each successive
secondary core around the primary protostar. The curve displays the instantaneous position of the secondary after completion of two orbital periods.
midplane temperature, and S1 2 = S 2 is the gas surface
plane is in turn partitioned into 150×150 square pixels, with10

ALMA 1.3mm観測
シミュレーション
The Astrophysical Journal, 857:40 (11pp), 2018 April 10

Sigalotti et al.

Figure 7. ALMA 1.3 mm continuum images of (a) disk A at 44.97 kyr, (b) disk B at 44.47 kyr, both at a viewing angle of 0◦, and (c) disk A at 44.97 kyr, (d) disk B at
44.47 kyr, both at a viewing angle of 45◦. The side of the disks nearest the observer is oriented toward the bottom of each frame. The source distance is assumed to be
at 100 pc from the observer. The angular resolution of each image is shown with the blue circle drawn in the lower left corner, corresponding to a synthetic beam sizeof 0″.05 (9 au at this distance). The simulation images were produced with an integration time of 1 hour. The color bar and numbers on the right indicate the ﬂux
intensity (in mJy beam−1).

at a distance of ≈15 au from the primary protostar at about
44.38 kyr, reaching a maximum separation of ≈19 au by
44.43 kyr. During this stage, the fragment accretes mass fromthe disk at a high rate as shown by the steep increase of MGI,
losing orbital angular momentum during this stage andmigrating to a minimum separation of ≈10 au by ≈44.6 kyrwhen it has already completed half of an orbital period.
Thereafter, as the disk dissipates, the accretion rate declines andthe secondary moves away in its orbit until it reaches amaximum separation of ≈18.2 au at ≈44.76 kyr by the time it
has completed a full revolution about the primary. The
secondary continues in its elliptic orbit and completes a second
period by ≈45 kyr when the calculation is terminated. At thistime, the mass of the secondary has become comparable to theprimary. During this second orbital period, the minimum andmaximum distances from the primary are ≈13.4 and ≈18.2 au,
respectively. This gives an orbital period of ≈225 years, whichcorresponds to a much shorter timescale compared to the Class0 lifetime of ∼30 kyr (André et al. 2000). The eccentricity was
found to decrease from ≈0.36 for the ﬁrst period to ≈0.151 forthe second period, suggesting that the secondary is approachinga more circular orbit as it comes close to the primary. However,sense, it is useful to simulate ALMA images to help interpretthe observations. Figure 7 shows simulated ALMA images ofdisks A and B at viewing angles of 0◦ and 45◦ and at two
different times during the evolution. These correspond to full
resolution images obtained by transforming the 2D emission
maps of Figure 6 with the aid of the simobserve andsimanalyze tools under Common Astronomy Software
Application for a wavelength of 1.3 mm (230 GHz; ALMA
band 6, conﬁguration C43-8 consisting of a full array of 50
antennas of size 12 m). A reference distance of 100 pc fromthe source was assumed. These images bear a remarkable
resemblance to the ALMA images of the Class 0 triple
protostar L1448 IRS3 (consisting of a tight binary and a more
widely separated tertiary) as observed when the tertiary, which
has the largest peak intensity in the system, is removed (Tobin
et al. 2016b). However, compared to this protostellar system,
the length scale in our results for disk B is smaller by a factorof ∼4.
3.4. Orbital Evolution of Second Protostar in Disk B
The orbital evolution of the second protostar around theprimary in disk B is shown in Figure 8. The secondary emerges
9

Figure 6. (a) Gas surface density (top row) and midplane temperature (bottom row) of disk A viewed face-on at four different times. (b) The same as (a) but for disk B.
By ≈44.97 kyr, when the calculation was terminated, disk A shows signs of fragmentation into a second protostar. The logarithm of the surface density Σ (in g cm−2)and the temperature T (in K) are plotted. Each row of images uses the color bars on their right sides. In all images, the central protostar is drawn with a black (surface
density) and dark blue (temperature) circle to increase the surface density and temperature contrasts.

11

18 March 7; in original form 2017 November 17

304 Dynamical evolu$on of stars and gas of young embedded stellar sub-clusters
Sills,1‹ Steven Rieder,2 Jennifer Scora,1 Jessica McCloskey1 and Sarah Jaﬀa3#52 Alison
MNRAS 477, 1903

ABSTRACT

We present simulations of the dynamical evolution of young embedded star clusters. Our initialconditions are directly derived from X-ray, infrared, and radio observations of local systems,
and our models evolve both gas and stars simultaneously. Our regions begin with both clusteredand extended distributions of stars, and a gas distribution that can include a filamentary
structure in addition to gas surrounding the stellar sub-clusters. We find that the regions
become spherical, monolithic, and smooth quite quickly, and that the dynamical evolution
is dominated by the gravitational interactions between the stars. In the absence of stellar
feedback, the gas moves gently out of the centre of our regions but does not have a significant
impact on the motions of the stars at the earliest stages of cluster formation. Our models atlater times are consistent with observations of similar regions in the local neighbourhood. We
conclude that the evolution of young protostar clusters is relatively insensitive to reasonablechoices of initial conditions. Models with more realism, such as an initial population of binary
and multiple stars and ongoing star formation, are the next step needed to confirm thesefindings.
Key words: galaxies: star clusters: general – galaxies: star formation.

さまざまな質量のbound clustersはどのようにして形成され進化するか？
Open clusters : liqle remaining gas, smooth and spherical stellar spaQal distribuQon
have had more than one episode of star formation (Gratton, Carretta
Younger embedded clusters: clustered or ﬁlamentary substructure& Bragaglia 2012).
この変化は、いつ、どのように起きるのか？
Gas成分の役割とその散逸の影響は？
f star formationis that ‘all stars
There are a number of as-yet
unanswered questions about the pro-te summary of the paper used tocesses through which embedded protoclusters evolve. For example,
da 2003) is実際に観測されたsub-structureをもつembeddedthat most Galactic star
most open clusters have little remaining
gas and dust,
andNGCare wellclusters（DR21,
M17,6357）を初期条
rotoclusters,
many of which will
described by a smooth, spherical stellar spatial distribution. Younger
件に採用して、星とガスの重力相互作用も入れて系の進化を追う
clouds as bound star clusters. Cerembedded clusters, on the other hand, more often show clustered12

Figure 4. Initial stellar positions for the DR21 simulations. Observed M
tIX stars are shown as red circles, and the inferred population of low-mstars are shown as small black points.

Figure 1. Initial conditions for DR21. X- and Y-axes correspond to RA Figure
and Figure
Figure3. Initialconditions
forand
NGC
6357.
See
Fig. 1toforRA
aThe
description.
The 3. Initial conditions for NGC 6357. See Fig. 1 for a description. The
1. Initialconditions
for
DR21.
XY-axescorrespondandlowest Figure
2. Initialconditionsfor M17.
See
Fig.
1 for
a description.
−2
Dec., respectively. Stars are shown as blueish-white objects; gas is included
Fig. 5levelis acorresponds
zoom-in of
mostdensity
centralsub-cluster,
lowestcorresponds
toa surfacedensity103
M⊙eachpc−2 , lowest
and contour
Dec., contour
respectively.
Stars
arelevelshown
blueish-white
objects;gas
included
to athesurfaceof 103M⊙ pc−2 , show
andlevelcontour
corresponds
to asa surfacedensityof 167
M⊙isofpc, and
0.5 higher.stars and their distribution
0.5 higher.
as diffuse emission, reflecting light from nearby bright stars. Gas is as
alsodiffuse
the sub-clusteringof thelow-mass
emission,fromnearbyeach subsequent
linelight
is 10
a factorof 100.5brighthigher.stars. Gas is alsoeach subsequentline is a factor
of 10
subsequent
linereflecting
is a factorofindicated with contour lines. The lowest contour level corresponds toaindicatedwith contour lines. The lowest contour level corresponds to aan ellipse determined by the observed stars. To show the ste−2 , and each subsequent line is a factor of
surface density of 150 M⊙ pc−2 , and each subsequent line is a factorsurfaceof planedensityof 150MTheof thesky).
are represented as objects with intensity
distributions in a slightly different way, Fig. 6 shows the cumula⊙ pcstars
1
100.5 higher. Figure created using FRESCO.1
100.5 higher.
Figurerelated
created using
FRESCO.
and colour
to theirmass and radius. Gas is rendered asradial distributions of the initial positions of the observed high-m

1908

A. Sills et al.

diffuse emission, and its smoothed surface density is also shownwith contour lines.
Table 2. Simulation parameters: stellar virial parameter (α vir ), velocity
dispersions
starspositions
and gas (σofstars
and inσ gas), fraction
In Figs
4 and of5,unclustered
we show the
stars
our
model ofof gas in each sub-cluster relative to stellar mass (fgas ), maximumDR21
velocity
thex–ysub-clusters
relative totoeach(Vrel ),planeand maximum
inoftheplane (equivalent
the other
RA–Dec.
of the sky).
offset of each sub-cluster from the z = 0 plane (Zoffset ).
Observed stars are shown as red circles, and the inferred low-massstars are shown as small black points.
Model nameα vir

Fiducial
Sub-virial
Super-virial
Low stellar sigma
High stellar sigma
Low gas sigma
High gas sigma
Low gas fraction
High gas fraction
Middle relative velocity
High relative velocity
Middle offset
High offset

1
0.08
2
1
1
1
1
1
1
1
1
1
1

σ stars (km s−1 )σ gas (km s−1 )fgas

Vrel (km s−1 )Zoffsetand inferred low-mass stars in DR21, measured in the plane ofsky from the centre of mass of our system. The radial distribution
the two populations are very similar until the edge of our simulat
area.
In Table 1, we give the global information about each of our melled regions (DR21, M17, and NGC 6357). These do not cha

MYStIX
Massive Young Stellar
complexes study in
Infrared and X-rays(Feigelson et al. 2013)

1
1
0
0
MNRAS1477, 1903–1912(2018)
from https://academic.oup.com/mnras/article-abstract/477/2/1903/4939295
1Downloaded1100
by National AstronomicalObservatory, Japan user
1on 25 May 2018
1
1
0
0
0
1
1
0
0
2
1
1
0
0
1
0
1
0
0
1
2
1
0
0
1
1
0.1
0
0
1
1
10
0
0
1
1
1
1
0
1
1
1
2
0
1
1
1
0
0.5 × zmax
11
10
1.0 × zmax
13
Figure 4. Initial stellar positions for the DR21 simulations. Observed MYSFigure 4. Initial stellar positions for the DR21 simulations. Observed MYStIX stars are shown as red circles, and the inferred population of low-mass
tIX stars are shown as red circles, and the inferred population of low-mass

describes the range of scales over which the fractal sub-structure
exists, with a low value indicating that the smallest sub-structures
are not much smaller than the overall cloud and a high value indicating that the smallest structures are very small compared to thecluster size. The density scaling exponent describes where most of
the mass is distributed in the hierarchy, with a high value indicatingthat most of the mass is confined to the smallest sub-structures anda low value indicating that the sub-clumps are not much more dense
than the larger clumps. Clusters that are not sub-structured but centrally condensed can be characterized by a radial density exponent
using the method described in Cartwright & Whitworth (2004) thatis incorporated into the Q+ analysis code.
What we found was that the fractal dimension of the DR21 fiducial simulation was consistently 1.58 (quite clumpy) until the subclusters merged at around 0.5 Myr. However, during that time, thenumber of levels decreases steadily and the density scaling exponent
Young sub-clusters
1909
increases, suggesting that the stars from the smallest sub-clump(s)are spreading out to fill the larger hierarchy. After the clumps merge
at ∼0.5 Myr, the region completely loses its fractal structure and isconsistent with a centrally concentrated system. The central stellar
density increases with time, as does the radial density gradient.
If we limit our analysis to only the high-mass stars (M >0.83 M⊙ ) that were taken from the MYStIX observations and exFigure 7. Snapshots of the first nine snapshots of the evolution of DR21
clude the low-mass stars added for completeness, we can compare
in our fiducial simulation, from top left to lower right. The snapshots aretaken 100 000 yr apart, with the first snapshot 100 000 yr after the simulation
the structure that would be theoretically seen by an observer to thebegins. Each snapshot is 5 pc on a side, otherwise as in Fig. 1. The subfull structure we have in the simulations. These results are shown inclusters have almost entirely merged
withinthe first million years ofclusterYoung
sub-clusters
1909
Fig. 8. For the initial conditions of DR21, this mass limited analysis
evolution, and the gas density has become quite spherical.
gave a lower fractal dimension and density scaling exponent thananalysis of the full cluster but a similar number of levels. This sugα vir = 2T/|W|. Our simulations have an additional gravitational
gests that the observable higher mass stars populate a more clumpy
potential from the gas mass, which effectively decreases the virial
distribution, but are not confined to the smallest sub-clumps. Theirparameter and shortens the time for collapse.
evolution showed the same trend as the full cluster. The high-massstars in M17 and NGC 6357 exhibit similar behaviour; the initialconditions have a low fractal dimension and high number of levels,
3.1 Sub-structure
but the sub-structure is quickly erased (within 1–2 Myr) and thecluster becomes centrally concentrated. It is interesting to note that
To characterize the sub-structure seen in our simulations, we used
when the full clusters are analysed, M17 would not be classified as+
the Q algorithm (Jaffa, Whitworth & Lomax 2017), which is not
fractally sub-structured and NGC 6357 only just so. It is only whenrelated to the virial parameter but rather characterizes the fractal diswe analyse the high-mass component of these two clusters that the+describes
sub-structured
the stars.
Q algorithm
guretribution
8. Figureofshowing
the Theevolutionof the cluster
structureaover
time.
Figure 9. Central stellar density versus core radius for all runs. The firstis detected. However, this difference in structure ise analyse
onlythemassive
stars
taken
from
the
MYStIX
observations,star cluster in terms of three parameters: the fractal dimension, the 7 Myr sub-structureof evolution was saved in snapshots every ∼100 000 yr, then +each
analysisambiguous
andmayas belimitationsin theThe
Q solid
this number
should beofsimilar
whattheobservers
be exponent.
able to detect.
if itcaused
were abyuniqueobservation.
levels,toand
densitywould
scaling
TheFilled
fractal di- density–radius pair plottedratherthananactualdifferenceinstructure
betweenthehigh-massmbols
indicate
lower
limits
(α
≥
3),
as
the
analysis
weperformcannotmension describes how clumpy or smooth the distribution is, with line is best fit to MYStIX observations. The various symbols and colours
and low-mass
imate
exponentsand
greaterthan 3.smooth.
The linesThebeginnumber
at the timeto differentcomponents.
simulations. Four of importance are the DR21 fidu1 radial
beingdensity
very clumpy
3 beingof levels correspond←10万年毎のスナップショット（5pc角の領域を示す）
By comparing our results to observaQons from the MYStIX
regions, we conclude that;
1. Our models match the observaQons well, and so theiniQal condiQons and choices of a priori unknown
parameters we used are realisQc and reasonable.
2. The subsequent evoluQon of these regions is quite
insensiQve to reasonable values of the parametersthat we varied.
We found that our iniQally very sub-structured systems
became spherical, monolithic, and smooth very quickly,
in both their spaQal and velocity distribuQons.
EvoluQon of embedded clusters at young ages appears to
be dominated by gravitaQonal interacQons between thestars, with only a small inﬂuence from the surrounding gas.

p in which the fractal structure is lost and the cluster becomes classified
centrally concentrated.

cial simulation (blue squares), the DR21 simulation with high gas fraction(orange squares), the M17 fiducial simulation (red triangles), and the NGC
6357 fiducial simulation (green circles). The grey squares are variations of
MNRAS
477,
1903–1912(2018)
ed from https://academic.oup.com/mnras/article-abstract/477/2/1903/4939295
the DR21 simulation with the parameters given in Table 2.
AstronomicalObservatory,
Japan user
2nal
Structuralproperties

y 2018e MYStIX observations provide a large sample of structural prop-14

018 March 7; in original form 2017 November 21

304 The diverse lives of massive protoplanets in self-gravita$ng discs
Stamatellos1‹ and Shu-ichiro Inutsuka2#53 Dimitris
MNRAS 477, 3110

ABSTRACTGas giant planets may form early on during the evolution of protostellar discs, while these are
relatively massive. We study how Jupiter-mass planet-seeds (termed protoplanets) evolve in
massive, but gravitationally stable (Q ! 1.5), discs using radiative hydrodynamic simulations.
We find that the protoplanet initially migrates inwards rapidly, until it opens up a gap in thedisc. Thereafter, it either continues to migrate inwards on a much longer time-scale or starts
migrating outwards. Outward migration occurs when the protoplanet resides within a gap withgravitationally unstable edges, as a high fraction of the accreted gas is high angular momentum
gas from outside the protoplanet’s orbit. The effect of radiative heating from the protoplanetis critical in determining the direction of the migration and the eccentricity of the protoplanet.
Gap opening is facilitated by efficient cooling that may not be captured by the commonly used
β-cooling approximation. The protoplanet initially accretes at a high rate (∼10−3 MJ yr−1 ),and its accretion luminosity could be a few tenths of the host star’s luminosity, making theprotoplanet easily observable (albeit only for a short time). Due to the high gas accretionrate, the protoplanet generally grows above the deuterium-burning mass-limit. Protoplanet
radiative feedback reduces its mass growth so that its final mass is near the brown dwarfplanet boundary. The fate of a young planet-seed is diverse and could vary from a gas giant
planet on a circular orbit at a few au from the central star to a brown dwarf on an eccentric,
wide orbit.
Key words: accretion, accretion discs – hydrodynamics – protoplanetary discs – browndwarfs – planetary systems.

last 20 yr have made possible thebe long when compared to the observed lifetimes of circumstellar discs (of the order of a few Myr; Haisch, Lada & Lada 2001;
Hernández et al. 2008; Muzerolle et al. 2010). Moreover, this model

15

Table 1. Simulation parameters.
Run idOpacity

a
αmin

b
αp,i
(au)

c
Mp,f
( MJ )

d
αp,f
(au)

efe

Run 1
Run 2
Run 3
Run 4
Run 5

BL94f
BL94f
BL94f
BL94f
SEM03g

0.1
0.1
0.2
0.1
0.1

40
40
40
40
40

24
15
29
25
31

40
17
37
36
53

0.035
0.0012
0.13
0.017
0.17

Notes. a Minimum SPH artificial viscosity. b Initial semimajor axis of the protoplanet. c Final protoplanet mass. d Final semimajor axis of the protoplanet.
e Final protoplanet eccentricity. f Bell & Lin (1994) opacity. g Semenov et al.
(2003)
3116opacity.D. Stamatellos and S.-i. Inutsuka

studies (Nayaksh& Montesinos 20
The way thatthe dynamical sta
difference betwee
protoplanet’s Hil
and the angular moutside the protop
mass of gas withthat moves towar

1 The same result is

the gas within dista

MNRAS 477, 3110–3135 (2018)

Downloaded from https://academic.oup.com/mnras/article-abstract/477/3/3110/4961145
by National Astronomical Observatory of Japan user
on 25 May 2018

Figure 3. The evolution of the mass Mp of the protoplanet for the runs
listed in Table 1. The increase in the protoplanet’s mass is initially rapid
as it opens up a gap and in most cases by the end of the simulation it has
become a brown dwarf (m > 11–16.3 MJ ; Spiegel, Burrows & Milsom 2011,
these limits are depicted by the horizontal dashed lines). When the radiative
feedback of the protoplanet is taken into account (Run 2), its mass growthis suppressed and its mass remains within the planetary-mass regime.

the gas that moves towards the protoplanet from outside its orbit(mout ) (Fig. 5, bottom). During the initial gap opening phase, in bothcases, gas is accreted on to the protoplanet mostly from the innerdisc. After the initial gap opening phase (∼5 kyr), we see that for
the run including the protoplanet’s radiative feedback (Run 2, red),
in which the disc is gravitationally stable because of the effect of the
Figure 2. The evolution of the disc surface density in the five simulations listed in Table 1. The star and the protoplanet as depicted as wide dots. Each rowprotoplanet’s feedback, most of the accreted gas comes from insidecorresponds to different snapshots of each run at times as stated on the graph. All runs show similar features, apart from the run in which the radiative feedbackthe protoplanet’s orbit. The total accreted gas in this case has lowerof the protoplanet is taken into account (Run 2; second row).

Figure 4. The ma
16 luminosi
accretionhigh, and thereforFigure 5. The difference between the mass of gas, within distance 1 to
1.2RH from the protoplanet, which moves towards the protoplanet frominside its orbit (min ) and from outside its orbit (mout ) (bottom) and the relateddifference between their angular momenta (Lin –Lout ). In the run including
the protoplanets radiative feedback (Run 2, red), most of the accreted masscomes from inside the protoplanet’s orbit, whereas in the case without
this feedback (Run 1, black) an almost equal amount of accreted masscomes from outside the protoplanet’s orbit. The average angular momentum
difference is positive for the radiative feedback run (for t > 5 kyr), but
negative for the standard (non-radiative feedback) run. After the gap opens
up in the disc, in the former case the protoplanet loses angular momentum
and migrate inwards, whereas in the latter case it gains angular momentum
and migrates outwards, while accreting higher angular momentum gas.

Figure 7. The semimajor axis evolution of an initially 1- MJ protoplanet
in a 0.1-M⊙ disc. The protoplanet initially migrates inwards on a Type I
migration time-scale (∼104 yr) but once a gap opens up, the inward migration stops and turns into an outward migration, except for the case when theradiative feedback of protoplanet is taken into account. In the latter case,the migration continues inwards but at a Type II time-scale (∼105 yr). (Thedetails of the runs are listed in Table 1.)

3118

D. Stamatellos and S.-i. Inutsuka

to accrete gas and eventually become a brown dwarf (unless the The estimated migra
disc dissipates by other processes, like e.g. photoevaporation, on in Table A1, and the mvalues are calculated aa faster time-scale). Nevertheless, when the protoplanet’s radiative The migration time-sc
feedback is taken into account the mass growth of the protoplanet protoplanet mass in Fig
has been suppressed and, at the end of the simulation, its mass is migration time-scales fby Ward (1997) and (Tlower by a factor of ∼2 than in the other runs.

et al. 2003).
The Type I migratio

4.2 Protoplanet migrationvI = −(2.7 + 1.1α% )The protoplanet initially interacts strongly with the disc resulting in
where α is the expon
fast inward migration, similar to the Type I migration for low-mass R −α% ), R% is the orbitalp4
planets in low-mass discs. The migration time-scale is only 10 yr of the disc at the positio
and it is similar for all five runs (see Figs 7, 8 and Tables A1, A2). position of the planet, aHowever, once a gap opens up the migration stalls, this happens at rate of migration of thedifferent times for each run. When the planet orbit plane is not the of the disc (e.g. Bate e
! "2
3
Hp
same with the disc mid-plane (Run 3), it is more difficult for theR
vI I = − α2Rpprotoplanet
up a gap.
In theat same
moreviscousFigureto8.open
The migration
time-scales
differentway,
times theduringthe runslisted isTable4)
1.theFilledboxesdifficult
correspond
empty
boxeset al. where α is the disc vi
the disc in
(Run
more
istotoinward
openmigration
up a gapand
(e.g.
Crida
to outward migration. The solid line corresponds to the Type I migration
actual Type II migratio
2006; Maliketal.
2015).
In
bothcases,
there
is
delay
in
openingtime-scale, whereas the dashed lines correspond to a combination of TypeI
up a gap
protoplanet
inwardsα (see(buttheonly imation (Crida & MorbandandType theII migration
assuming migrates
different discfurther
viscosity parameter
Kley 2015, 2017). For
textcomparedfor details). Oncegap otheris openedup, interactionsof the runs
protoplanetby ∼4 au
withthetheruns).
In the four
without that due to gas crossin
with the gravitationally unstable gap edges lead to outward migration. Whenradiativethefeedback
from thetheintoprotoplanet
switchesradiative feedback
fromprotoplanet,
the planet is taken
account the gap
edges to faster or slower. The tomigration.
Ataretheend
simulations(20 kyr), and Type II regimes, isare hotter, andtherefore theystableandofthetheinwardmigration continues
Figure 18. The evolution of a protoplanet in a 0.1-M⊙ protostellar disc. The protoplanet is placed at different initial orbital radii within the disc: 5, 10, 20, a slow outward
but at a longersemimajor
time-scale. axis has moved to 40, 38, and 36 au,
40, and 80 au (top to bottom row; sRun1, . . . , sRun5, respectively). There is a difference in the migration pattern depending on the position of the protoplanet the protoplanet’s
17 vI
+
v=
in the disc.
1 + (Mp /Mt )3for Runs 1, 2, and 3, respectively, from an initial semimajor axis of

3048; inTheforma$onof rings
18 March
original
form 2017 December
17 and gaps in magne$cally coupled disc-wind systems:#54 ambipolar diﬀusion and reconnec$on

Scoq S. Suriano,1‹ Zhi-Yun Li,1 Ruben Krasnopolsky2 and Hsien Shang2
MNRAS 477, 1239

ABSTRACT

Radial substructures in circumstellar discs are now routinely observed by Atacama Large
Millimeter/submillimeter Array. There is also growing evidence that disc winds drive accretionin such discs. We show through 2D (axisymmetric) simulations that rings and gaps develop
naturally in magnetically coupled disc-wind systems on the scale of tens of au, where ambipolar
diffusion (AD) is the dominant non-ideal magnetohydrodynamic effect. In simulations where
the magnetic field and matter are moderately coupled, the disc remains relatively laminar withthe radial electric current steepened by AD into a thin layer near the mid-plane. The toroidal
magnetic field sharply reverses polarity in this layer, generating a large magnetic torque thatdrives fast accretion, which drags the poloidal field into a highly pinched radial configuration.
The reconnection of this pinched field creates magnetic loops where the net poloidal magnetic
flux (and thus the accretion rate) is reduced, yielding dense rings. Neighbouring regions withstronger poloidal magnetic fields accrete faster, forming gaps. In better magnetically coupledsimulations, the so-called avalanche accretion streams develop continuously near the discsurface, rendering the disc-wind system more chaotic. Nevertheless, prominent rings and gaps
are still produced, at least in part, by reconnection, which again enables the segregation of
the poloidal field and the disc material similar to the more diffusive discs. However, thereconnection is now driven by the non-linear growth of magnetorotational instability channel
flows. The formation of rings and gaps in rapidly accreting yet laminar discs has interesting
implications for dust settling and trapping, grain growth, and planet formation.
Key words: accretion, accretion discs – MHD – protoplanetary discs – ISM: jets and outflows.

18oh0.26
oh2.6
oh26
beta1e2
beta1e4

0.922
0.922
0.922
0.0922
9.22

4.408
4.408
4.408
4.408
4.408

0.25
0.25
0.25
0.25
0.25

9.71
9.71
9.71
97.1
0.971

2.5
25
250
–
–

Figure 1. A representative (‘reference’) axisymmetric simulation. Shown are the mass volume density (logarithmically spaced colour contours in units ofg cm−3 ), the poloidal magnetic field lines (magenta), and the poloidal velocity unit vectors (black). Panels (a)–(d) corresponding to simulation times of 0, 200,
1000, and 2500 inner orbital periods, respectively. (See the supplementary material in the online journal for an animated version of this figure.)

3 REFERENCE MODEL AND THE
F O R M AT I O N O F R I N G S A N D G A P S T H RO U G HRECONNECTION(left to right), where t0 = 1 yr is the orbital period at the inner edgeof the simulation domain (r0 = 1 au). The disc-wind launching
19 from
proceeds in an inside-out fashion, i.e. the wind is launchedlarger disc radii as the simulation progresses. By 1000t0 (one orbit

Figure 3. The reference simulation at time t/t0 = 2500. The colour contours show (a) the logarithm of poloidal velocity (cm s−1 ), (b) the logarithm ofplasma-β, (c) the ratio of the toroidal to the poloidal magnetic field components, Bφ /Bp , (d) the logarithm of density (g cm−3 ), (e) an axisymmetric, face-on
view of the disc surface density distribution normalized to the initial power-law distribution, % i = % 0 (r/r0 )−1/2 , and (f) the mass flux per unit polar angle,
dṀ/dθ = 2πr 2 ρvr sin θ , normalized to Ṁ0 = r02 ρ0 cs,0 . Poloidal magnetic field lines (i.e. magnetic flux contours) are shown in grey in panel (c). Panel (d)shows two specific poloidal magnetic field lines with mid-plane footpoints at r = 8 au (magenta; see Fig. 4) and r = 7 au (black; see Fig. 5). Poloidal velocity
unit vectors are plotted in black in panels (a), (c), (d), and (f). (See the supplementary material in the online journal for an animated version of this figure.)

20

1248

S. S. Suriano etRing
al. and gap formation

1247

compared to its poloidal component. The now weaker poloidal fieldstrength can be seen clearly in the last panel from the lack of fieldlines in the post-reconnection region.
3.3.3 Reconnection-driven ring and gap formation
The general picture for the reconnection-driven ring and gap formation is as follows. A poloidal magnetic field line that initiallythreads the disc rather smoothly is dragged by preferential accretion near the mid-plane into a highly pinched radial configuration(see the upper-right panel of Fig. 7 for an illustration). Reconnection of the highly pinched field line produces a poloidal magnetic loop next to a poloidal field line that still threads the discmore smoothly (see the lower-left panel of Fig. 7). After reconnection, the material trapped on the magnetic loop is detached fromthe original (pre-reconnection) poloidal field line, which resultsin the separation of matter and (poloidal) magnetic flux. SpecifFigure 7. Poloidal magnetic field lines at four different times are shown in grey, with a reconnecting field line highlighted in black. The colours show the ratio
ically, there is no net poloidal magnetic flux passing through theof the toroidal to poloidal magnetic field, Bφ /Bp .
matter on the loop, and there is less matter remaining on the original poloidal field line (since some of the matter on the original field
diffusivity, such as resistivity), prominent avalanche streams would
3.3.2 Reconnection of pinched radial magnetic field
is now on the detached magnetic loop). The net effect is a redishave developed near the disc surface, which would driveline
the entire
The secondary component of the mid-plane current sheet develops
disc and its envelope into an unsteady state and make tribution
the forma- of the poloidal magnetic flux away from the loop-forming
as a result of mass accretion in the disc, which is concentrated intion of a laminar mid-plane current sheet impossible (seePaperintoI
regionanadjacent
regionwhere
the poloidal
flux accumulates.
near the mid-plane(see Fig. 6).
theprimary radialcurrent(Jr ) sheet
and Section 4.1 below).1 Secondly, as first stressed by BrandenThe mass accretionisdrivenby the removalof angularmomentumSince
mass
accretion
tends
to
be
fasterinregionswithstrongerburg & Zweibel (1994), AD tends to steepen the magnetic gradientwhich
due to the Lorentz force (∝Jr Bθ ) in the azimuthal
poloidalthat thedirection,
reconnectionnear a magnetic null, i.e. where the magnetic field changes
polarity. magnetic fields, it is not surprisingpeaks in the radial current sheet where the toroidal magnetic field
The reason is that the Lorentz force associated with theinduced
magnetic variation of the poloidal field strength with radius wouldchanges sign. Pictorially, as the φ component of the magnetic fieldpressure gradient drives the ions (relative to the neutrals) towards
lead to a spatial
variation
massbelow
accretion
would inchangesfrom thein+φthedirection
the discrate
to thethat−φ direction
the null from both sides. Since the field lines (of opposite polarabove thedisc in a thin
mid-plane
the magnetic
fieldi.e.
linestheturnlead to spatial
variationin themass layer,(surface)
density,
ity across the null) are tied to the ions, they are dragged
towards
become severely kinked
in the azimuthal direction. The sharp field3
the null as well, leading to a steepening of the magnetic
gradiformation of ringsand gaps.
kink generates a magnetic tension force in the −φ direction thatent, which in turn yields a stronger magnetic force that drives the
The preferential
poloidallines
insideexerts a concentration
braking torque andof
drivesthe discfieldaccretionin thecurrenttheions and the field lines even closer to the null. Since the ambipolar
sheet. of the disc is illustrated pictorially in Fig. 8(a),low-density
gaps
magnetic diffusivity, ηA , is proportional to the field strength (seeThe accretion through the mid-plane current sheet now drags theequation 5) and thus vanishes at the null, this steepening
would
which
is a snapshotof the inner part of the reference simulation
poloidal magnetic field lines inwards with the flow. This results inresult in an infinitely thin, singular, current sheet in principle.
In(up to 15 au)a pronounced
at a representative
time
when magnetic
the simulationhasradial pinch of
the poloidal
field, which
practice, the thickness of the current sheet is limited by the gridtransports thestate,poloidalflux The
inwards
and concentration
yields anotherfieldreached a quasi-steadyt/tmagnetic
0 = 2500.
resolution.
in the azimuthal direction, Jφ . Eventually, theis quantified
incurrent
panelcomponent(b), whichplots the vertical component of the
Sharp magnetic field reversals that give rise to thincurrent
radial pinch becomes so severe that the magnetic field reconnects,
sheets are prone to reconnection. However, this is not thecaseforred) atthe mid-plane
andaresurface
densitymagnetic fieldforming(Bz , inpoloidal
magnetic
field loops that
reminiscentof the ofthe mid-plane current sheet shown in Fig. 6, because it is protearing bothmode relative
instabilityto(Furthal. 1963;
see also&
the disc (in black),
theiretinitialvaluesat tZweibel= 0. Two
duced by the reversal of the toroidal field component Bφ , andYamada 2009). An example of the reconnection process is shown
featuresit is impossible to reconnect oppositely directing toroidal
fields stand out: (1) there is a strong anti-correlation between
in Fig. 7, which plots lines of constant poloidal magnetic flux, i.e.
curunder the adopted axisymmetry.2 Nevertheless, this primary
the surface
densityverticalfieldstrength,
as expected
in tothepoloidal andmagnetic
field lines,
alongwith the ratioof the toroidal
rent sheet leads to another current component that does allow
reconnection-induced
scenario
ring
and
gap
formation,
and(2)
contours),
at a radius centred
poloidal magnetic
field, ofBφ /Bp (colour
for reconnection.
230. At
the first frame
shown,
on r = 6.5 aufieldfrom time
t/t0 = 200 to
the poloidal magnetic
is typicallymuchweakerin thehighthe poloidal magnetic field highlighted in bold already has a kinkclose
to zero) than in the lowdensity ring regions
z ispinch
across the (where
mid-plane.BThegrows with time until the field lineisincreasedoveritsapproximately
initial value
densitygap
regions(whereBisstretchedalmostparalleltothemid-planeoverfourbyz
1 In this reference run, AD does not appear to suppress the development ofradial grid cells at t/t = 210. By the next frame shown at t/t0 = 220,
factorflux
the MRI completely. Channel-like features are evident at largearadii
whereof 2–4). This drastic 0segregation of poloidal magnetic
the field has reconnected, forming a loop near r = 6.4 au. In thethe disc is better coupled to the magnetic field as measured by relative
the Elsasserto matter
appears
to
behas
quite
stable withtimethe(persistingforlastframe,
theloop
disappeared;
however,
processof
its
number (see Fig. 1c).
of
innerorbits),despite
the
fact
that
there
is
continued
2 Reconnection of the highly pinched toroidal field is expectedthousands
formation
has
left
a
lastingmarkonthemagnetic
field
structure.
in 3D, andTheregion now
has athemuchlargerand
toroidalmagnetic field component
will be explored in a future investigation.
mass accretion
through
both
rings
gaps.

Figure 8. The density, magnetic field strength, and velocities plotted up tor = 15 au at t/t0 = 2500. The panels show (a) the logarithm of the density(g cm−3 ) in colour and the poloidal magnetic field lines (i.e. magnetic flux
contours) in black, (b) the surface density normalized to the initial radial
distribution, ! i ∝ r−1/2 , and the vertical magnetic field strength at the
mid-plane normalized to its initial distribution, Bz, i ∝ r−5/4 , (c) the radial
velocity (km s−1 ) of neutrals (black) and ions (red) at the mid-plane, and(d) the density-weighted vertical average of the radial velocity (km s−1 )of neutrals (black) and ions (red) over z = ±2h. (See the supplementarymaterial in the online journal for an animated version of this figure.)

presence of significant AD, this is not necessarily the case. Indeed,at the time shown in Fig. 8, the ions are moving outwards in severalof the gaps, especially the two near r = 5 and 6.5 au. The ionsare forced to expand against the infalling neutrals by an outward
Lorentz force due to either a temporary poloidal flux concentrationnear the inner edge of the gap [in a manner that is reminiscent
21ofthe forced ion–neutral separation in the AD shock in magnetizedReceived September ....; accepted ...

304 Unveiling the remarkable photodissocia$on region of M8
Tiwari1,⋆,, K. M. Menten1, F. Wyrowski1, J.P. Pérez-Beaupuits1, 2, H. Wiesemeyer1, R. Güsten1, B. Klein1, C. Henkel1, 3#55 M.
A&Ap in press
ABSTRACT

Aims. Messier 8 is one of the brightest HII regions in the sky. We collected an extensive dataset comprising multiple submillimeter spectral lines from neutral and ionized carbon and from CO. Based on it, we aim to understand M8’s morphology and that ofits associated photo dissociation region and to carry out a quantitative analysis of the regions’ physical conditions such as kinetic
temperatures and volume densities.
Methods. Using the Stratospheric Observatory For Infrared Astronomy (SOFIA), the Atacama Pathfinder Experiment (APEX) 12 m
and the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescopes, we have performed a comprehensive imaging survey ofthe emission from the fine structure lines of [C II] and [C I] and multiple rotational transitions of carbon monoxide (CO) isotopologues
within 1.3 ⇥ 1.3 pc around the dominant Herschel 36 (Her 36) system composed of at least three massive stars. To further explore themorphology of the region we compare archival infrared, optical and radio images of the nebula with our newly obtained fine structureline and CO data, in particular with the velocity information they provide. We perform a quantitative analysis, using both LTE andnon-LTE methods to determine the abundances of some of the observed species as well as kinetic temperatures and volume densities.
Results. Bright CO, [C II] and [C I] emission has been found towards the HII region and the photodissociation region (PDR) in M8.
Our analysis places the bulk of the molecular material in the background of the nebulosity illuminated by the bright stellar systems
Her 36 and 9 Sagitarii. Since the emission from all observed atomic and molecular tracers peaks at or close to the position of Her
36, we conclude that the star is still physically close to its natal dense cloud core and heats it. A veil of warm gas moves away fromHer 36 toward the Sun and its associated dust contributes to the foreground extinction in the region. One of the most prominent star
forming regions in M8, the Hourglass Nebula, is particularly bright due to cracks in this veil close to Her 36. By using radiativetransfer modeling of di↵erent transitions of CO isotopologues, we obtain H2 densities ranging from ⇠ 104 – 106 cm-3 and kinetic
temperatures of 100 – 150 K in the bright PDR caused by Her 36.
Key words. galactic: ISM — galactic: individual: M8 — submillimeter: ISM — ISM: HII region — ISM: clouds

Introduction

e influence of bright stars on their surrounding interstellar
dium is immense. Their strong ultraviolet and far-ultraviolet
lds give rise to bright HII regions and photodissociation reons (PDRs). These are the best grounds to study the e↵ect of
V and far-UV photons on the heating and chemistry of the instellar medium (ISM). The fine structure lines of C+ and O,

ous molecular and atomic species allows us to derive the physi
conditions in PDRs.
Messier 8 (M8) is located in the Sagittarius-Carina arm, n
our line of sight towards the Galactic Center. It is located at a dtance ⇠1.25 kpc (10 corresponds to 0.36 pc) from the Sun (Damani et al. (2004); Arias et al. (2006)) with an error of ⇠ 0.1 k
22(Tothill et al. 2008) and is about 34 ⇥ 12 pc in diameter.
M
is associated with the open young stellar cluster NGC6530,

9 SgrHer 36

1.5 pcFig. 1. The left panel shows the Lagoon Nebula, Messier object 8 (M8) or NGC 6523, in the constellation of Sagittarius, as seen by the Kitt Peak
4-meter Mayall telescope in 1973. North is up, east to left. Credits to "National Optical Astronomy Observatory/Association of Universities for
Research in Astronomy/National Science Foundation" Copyright WIYN Consortium, Inc., all rights reserved. The right panel shows the WISE
3.4 µm image of the M8 region around Her 36 (marked with a star) investigated in this paper. The contour levels are 10% to 100% in steps of 10%of the peak emission ⇠ 4000 data number. The 1.5 pc indicated in the lower left correspond to 40 .

Fig. 4. Colour maps of the integrated intensity (left to right) of the J = 1 ! 0 transition of 12 CO, 13 CO and average of H↵40, 41,42 and 43 lines
toward Her 36 which is the central position ( ↵ = 0,
= 0) at R.A.(J2000) = 18h 03m 40.3s and Dec.(J2000) = 24 220 4300 , marked with antralbaselines ofThefirst orderand subsequently
produceddata cubes
wassteps of 10% of the corresponding peak emission given in Table 1.
to the
SOFIA
observations.
The
pointing(< 300 ) inasterisk.
contourlevels
are
10%
(>
3
⇥
rms,
given
in
Table
1)accuracy
to 100%with beam/3 sampling. For the beam sizes, see Table 1.
maintained by pointing at bright sources such as RAFGL5254
R Dor shown
every 1 – 1.5
forward efficiency
Fe↵ each= 0.95 map.
All maps are plotted using original beamandsizes
in hrs.
theA bottom
left of2.2. The APEX data

was used for all receivers, and the beam coupling efficiencies
Be↵ = 0.62, 0.69, 0.63, 0.43 and 0.32 were used for the PI230,
FLASH+ 340, FLASH+ 460, CHAMP+ 660 and CHAMP+ 810 receivers, respectively.

Table 1. Line parameters of observed transitions.
12

13

Observations of low- and mid-J CO and CO transitions were
performed with the Atacama Pathfinder Experiment (APEX)
12m submillimeter telescope (Güsten et al. 2006) during 2015
June - August and 2016 July, September and October. As shownin Table 1, we used the following receivers: PI230 to map the
low-J 12 CO and 13 CO transitions, FLASH+ in the 345 and 460
GHz bands to map the mid-J CO transitions and CHAMP+ in
low and high frequency sub arrays to map the higher frequency
mid-J CO transitions.
We used the PI230 receiver to map C18 O J = 2 ! 1 at
219.560 GHz , 13 CO J = 2 ! 1 at 220.398 GHz and 12 CO J
= 2 ! 1 at 230.538 GHz. FLASH+ was used in the 345 GHzband to map the 12 CO J = 3 ! 2 transition at 345.795 GHz.
FLASH+ was also used in the 460 GHz band to map the 13 CO J
= 4 ! 3 transition at 440.765 GHz, C18 O J = 4 ! 3 at 439.088
GHz and [C I] 3 P1 ! 3 P0 fine structure line at 492.160 GHz. The
CHAMP+ receiver was used to map 12 CO J = 6 ! 5 at 691.473
GHz and 13 CO J = 6 ! 5 at 661.067 GHz in the low frequency
sub array complemented by 12 CO J = 7 ! 6 at 806.651 GHz
and 13 CO J = 8 ! 7 at 881.272 GHz in the high frequency sub
array.
All maps shown in Fig. 3 were observed in on-the-fly(OTF) total power mode centered on R.A. 18h 03m 40s .3; Dec.
24 22’4300 (J2000), which corresponds to the position of Her
36. The maps obtained from the observations done in July,
September and October 2016 have a size of 24000 x 24000 . Mapsthat were obtained from the observations done in 2015 are comparatively smaller in size. We integrated 0.7 s per dump for allmaps and enough coverages were performed to reach the rms
noise levels as mentioned in Table 1. The o↵set position relative to the centre at (300 , 300 ) was chosen as reference, similar

Transition

J=1!0
J=2!1
J=3!2
J=6!5
J=7!6
J = 11 ! 10
J = 13 ! 12
J = 16 ! 15

2

rms (K km s 1 )TelescopeCO
115.271
0.73
22.5
610.2
230.538
0.65
28.7
355.5
345.796
0.73
19.2
210.2
691.473
0.43
9.6
580.1
806.652 All maps0.34
673.4
shown in Fig. 48.2
were observed in OTF total power
centered on Her 36. Each subscan lasted 25 s and the in1267.014 mode
0.68
22.9
130.9
tegration time on the o↵-source reference position was 5 s. Theo↵set position
relative to the19.1
centre at (30 , 30 ) was similar 155.3
to1496.923 that
0.68
used for the SOFIA and APEX mapping and the pointing3 ) was maintained
1841.346 accuracy (< 0.70
14.8by pointing at the bright cali-46.2

1.1
0.4
0.8
4.0
19.0
3.5
2.4
1.8

IRAM 30m/EMIR
APEX/PI230
APEX/FLASH+APEX/CHAMP+
APEX/CHAMP+SOFIA/GREATSOFIA/GREATSOFIA/GREAT

CO
0.73
23.5
64.9
0.65
30.1
92.6
www.iram.es/IRAMES/mainWIKI/IRAM30mEfficiencies
0.59
15.0 Article number, page 3 of 198.1
18
0.45
10.0
158.8

0.4
0.6
2.5
3.2

IRAM 30m/EMIR
APEX/PI230
APEX/FLASH+APEX/CHAMP+

12.7

0.6

APEX/PI230

34.0

1.8

APEX/FLASH+

728.5

2.1

SOFIA/GREAT

11.2

0.2

IRAM 30m/EMIR

J=1!0
J=2!1
J=4!3
J=6!5

110.201220.399440.765661.067

J=2!1

219.5613

492.160

P1 ! 3 P022.3. The IRAM 30 m data

Frequency (GHz)⌘mb

Peak line flux (K km s 1 )✓mb (00 )Observations of low-J 12 CO, 13 CO and hydrogen recombination
12 with the IRAM 30 m telescopeline observations were performedin August 2016. We observed the whole 3 mm range using the
EMIR receivers (Carter et al. 2012). We simultaneously mappeda region of 24000 x 24000 , which is similar to the size of most othermaps previously observed with SOFIA and APEX, in the 12 CO
and 13 CO J = 1 ! 0 transitions (Table 1) and the hydrogen recombination lines H40↵ to H43↵ at 99.023 GHz, 92.034 GHz,
85.688 GHz and 79.912 GHz, respectively. Molecular high density tracers, also detected in our wide spectral band observation,
will be analyzed in a subsequent paper.

0

0

00

brator 1757-240 every 1 – 1.5 hrs. A forward efficiency Fe↵ =0.95 and a beam coupling efficiency
Be↵ = 0.69 were adopted
13
for the EMIR receivers. These values were taken from the latest(2015) commissioning report4 .

4

0.65

C18 O
30.2
12

0.59

12

P3/2 ! P1/2

1900.53

0.70

H40 – 43↵

80–90

0.73

C

13.5
C

+

14.8
H↵
⇡ 26

23

M. Tiwari et al.: Unveiling the remarkable photodissociation region of M8

Fig. 8. Line profiles at di↵erent o↵sets (00 ) relative to Her 36 are shown in di↵erent cdetailed
positions,intensity
see Sect.of:3.4.
panels:!II5] 12
COthe
andJ =
J =6 1!!5 0
Fig. 5. Scatter plots and correlation coefficients r between thethe
velocity
integrated(a) Upper[C II] and[C IJ] =(b)6 [Candof [Cdata
I] and[C IIwere].̃ Allextracted
spectra were
their originalbeam
siztransition of 12 CO and (c) [C I] and the J = 2 ! 1 transitiontransitionsof 13 CO. All
pointsfromextracted
velocityfrom
integrated
intensity
maps
convolved to the same beam size of 3100 .

M. Tiwari et al.: Unveiling the remarkable photodissociation region of M8µm (bandcorresponds to the north-east of Her 36 where [C II] is morecloserexto
12
tended. The lower right, where CO J = 6 ! 5 emission
inten13 km s-1
originatessifies at a faster rate than [C II], corresponds to the south-westof sphologyO↵set (00 )a V (km s 1 )V (km s 1 ) T peak b(K)Her 36, where 12 CO J = 6 ! 5 is much more prominent.
The
HII
regio
12
CO J = 6 ! 5
correlation of [C II] with [C I] has a correlation coefficient of r
( 40, 35) 10.33 (0.06)
3.62(0.16)
77.24
s-1
10 km
= 0.473 and again shows two di↵erent branches corresponding
warm 3.30 (0.26)
5.92 (0.10)
58.41
to di↵erent regions. The upper left, where [C II] emission
3.4. getsSpe( 13, 8) Sgr 9
PDR
10.43 (0.04)
2.92 (0.10)
134.43brighter for an almost constant [C I] emission, correspondsdiffeto
12
13
HII region
6.39(0.03)
2.42 (0.07)
45.67
theobtained
north-eastfromof Her
36. The modelingbranch in theFig. 13. Results
RADEX
forlower
COright,
and similar
CO atto t(0, 0) 2 km s-1
00 5 (b), corresponds
10.42 (0.02)
3.58 cold(0.04)
97.94
Fig.
8 sho
Fig.
the south-west
and dense
000 ) and at thethe[C situation
II] peak (shown↵ = in30.0,= 2.000 ). Ourtoobserveddatap
Her 36(30, 2) 10.28 (0.09)
4.73 molecular(0.20) cloud
62.60[CcorreI] andFig. 2. Colour maps of the integrated intensity of the [C II] 158 µmand J = 11 ! 10, J = 13 ! 12 and J = 16 ! 15 transitions of CO toward00
is wellof Her
36.same
In contrast
to thesecorrelations,[C I]densities
maps
convolved
totheresolution
of
31
.
The
columnused
Her 36 which is the central position ( ↵ = 0,
= 0) at R.A.(J2000) = 18 03 40.3 and Dec.(J2000) = 24 22 43 , marked with an asterisk. The
13(60, 27)
11.49 (0.07)
2.51 (0.17)
53.77
36. Line
latedresultswith obtained
CO J = 2by!varying
1 with rthe
= 0.908.
This
resembles
thecontour levels are 10% (> 3 ⇥ rms, given in Table 1) to 100% in steps of 10% of the corresponding peak emission given in Table 1. All maps are
7 km s-1(a) and (b) show
kinetic
temperaturefromplotted using original beam sizes shown in the bottom left of each map.
13
13
case
M17towards
SW, fora[Ccorrelation[CinI]Table
with
Fig. (d)9. LineprofilesHerwhich
36byforvarying
I] 3 P1 ! 3coefficient
P0 ,density
CO J of=
4
CO J = 1 ! 0(c) andshow
results
obtained
the
H
from
10
–1
2
6! 5, 12 CO13JCO
= 6!
5 and
] 2 P3/2! 2 P1/2to. All
are extracted velocity
J=
2 ![C1IIwas
reportedbe spectra
0.942 (Pérez-Beaupuits
et al.
ObserverAll data reduction was done using the CLASS and MIRA ous transitions of 12 CO, 13 CO and C18 O are presented in Ta( 40, 35) 8.47 (0.04)
3.00 (0.11)
8.26
fixedfromat 120
K.
00
J
=
6
!
5
12
maps
which
were
convolved
to
the
same
beam
size
of
15
.
A&A
proofs:manuscriptno.
main
2015c).
programs that are a part of the GILDAS software package and ble 1. Fig. 2 shows velocity integrated intensity maps of the CO( 13,8) Her8.76(0.03)
all the observations are summarized in Table= 11 ! 10,intensities
13 ! 12 andnormalized
16 ! 15 transitions.
emissionFig. 1.15. Left: VelocityJintegrated
to theirThevalue
at peak vs o↵set(00 ) from36 along
the green2.58
arrows(0.08)
shown in the12.61
right.
of the ge
00 peaksin allthe
the secondlines hasATLASGAL
a similar spatialpeak
distribution
and
2.68
15.16
They follow a path fromat ( 53, 2300 )areto the second(0,[C 0)II] peak at8.91(6000(0.01), 2700 ) via Her
36 at(0.03)(000 , 000 ). Plotsarediscussed
found at about 12
the same o↵set position18( ↵ = 5.000 ,= 5.000 )3. thePhysicalparameters
calculated image
from rotational
shown for [C II], [C I],north-west
J = 6 !of5Her
CO,J = 2 ! 1 C O, H↵, ATLASGAL
870 µm(30,and WISE
emission;
Schematic
diagram 11.15
ofTablethe with
3. Results2) 3.4 µm
9.37(0.03)Right:2.48(0.09)ATLASGAL
dust continuum
that peaks diagrams.
at Her a curved
36.
prominentofM8
pertinenttothediscussioninthispaper.
The
cold
and
dense
molecular
cloud
is
in
the
background
shown
in
blue.
3.1. Peak intensities of the molecular line
emission optical features
3.2.
Channelmaps
Figs. 3 and 4 show velocity integrated intensity maps of low(60, 27)
10.48 (0.02)
1.45 (0.05)
8.46 36 and also traces the cold molecular cloud in the north-west. east of H
13
The foreground gas ofandthemid-J
warm
PDR veilrecedingaway
Her(⇠ 9 km s 1 ) with lower velocities (2 – 7 km s 1 ), while the HII region istransitionsof 12isCO,
CO andC18 O,from
i.e. theJ =36
1!
The maxima of the distributions of the velocity integrated inThe dust emission morphology
is similar to the 12 CO, 13 CO and correspon
12 3
13 CO
powered by both stellar
9 Sgr.
0, 2systems,! 1, 3 !Her2, 636! and
5, 7 !
6 transitions of 12 CO, the J = 1
C P1 ! 3 P0
Optically thin
Correction
factorincludedtensities of the emission in the [C I] and [C II] lines and vari18
13! 0, 2 ! 1, 4 ! 3 and 6 ! 5 transitions of CO and the J = 1
CO
distribution.
Secondly,
weusedthedatafrom theNational
Ra- at ( ↵of=In order 13to investigatedi↵erencesinthe distribution! 0 transition of C18 O. The intensities of the low-J transitions
www.iram.fr/IRAMFR/GILDAS/( 40, 35) 9.40 (0.63)
5.91 (2.16)
2.84 dio Astronomy-2 )13 CO)-2a(K)a(K)
Observatory(NRAO)/
Very
Large
Array
(VLA)
T
N(
CO)
(cm
T
N(
(cm
2
2 of trot ionized and atomic carbon, channel
rot
peak
maps
of
the
P3/2) !
P1/2
18)s . Therefore
9.47 (0.25)
3.63that(0.63)
5.62 Archive Survey (NVAS6 ). Fig.
4.4. CO, [C I] and [C II] luminosities
veil seen( at13,
7 km
we suggest
the low velocity
Article number, page 4 of 18
7
(b)
shows
the
[C
II] velocity 16
in3
3
16transition3.6
of [CII] are comparedto those of[C I] P1 !=P02.
8.3 ± 1.3⇥ 10
8.6 ± 1.5
6.2the⇥ 10[C II] probes
directlythegas(0.18)of the veil 3.86
that forms
a foreground
(0,
0)
9.89
(0.47)
7.65
intensity
mapInoverlaid
with thecm radio continuum
We obtained the total luminosities of the CO spectral line energy PDR illuminated by Her 36. On a fainter level, weak emission tegrated
mid-J
transition.
be1.3
emission[C IICO
]417 ± 3.7
3.6 ⇥Fig.
10166, it can 30±seen
4.8 that the5.4⇥ 1016 from
(30,
2)
10.05
(0.30)
4.01(0.64)
5.45
12
peaks
veryisclosetoHer 36
andcomparedtraces free-free
distributions (SLED), and of the [C I] 609 µm and [C II] lines, from this veil is also seen in the CO maps at low velocities (5 image whichCO
16
16 from(lowerpanels)spread
that
[Cand
I]
90 ± 21.6
1.9 ⇥
10more39 ±out
4.9asThis
4.6 ⇥to10HII
11.80 (0.09)
0.68 (0.35)
5.85 emission
over the total observed region seen in the maps of Figs. 2, 3 and – 7 km (60,
from
the
HII
regionNGC6523/33.
compact
1 27)s ). This foreground veil is receding away from Her(upper panels). In the velocity
2 km
s 161 to 6ATLASG
km s 1±range10 from2.1⇥ 10
12 + 24, as derived by Solomon et al. (1997); Carilli & Walter (2013) 36 towards us and to the west
bykm
thes H184
recombinationlines
with [CtheC withP3/2! 2 P1/2
1
a changein the line-of-sight region, which
west tand is15alsokm straced
to 17
there
is no emission
fromI] whileand scaled it for galactic sources:
IRAM there
30m istelescope,also[Cshown
in Fig. 4 (right panel).
⌧ isfrom
Transition⌧
velocity. This is consistentboth high5.40 with(0.14)
2.84velocity(0.30) red shifted
33.45the Transition
II] close to Her 36 and towards the east
emission40, 35)blue shifted emission of H↵, [N II] and [S II] Thirdly, the Wide-field Infrared Survey Explorer (WISE) imaged
and low( velocity(4)
L = 1.04 ⇥ 10 9 S V⌫D2L ,9.76 (0.11)
3.80 (0.31)
45.05 J = 1 !of0 it, respectively.
1.68 In the Jrange= 1 !from
0 7 to 9 1.92
km s 1 both [C]A IIlowdoublets, [O III], and absorption lines of the sodium D doublet the sky at four mid-infrared wavelengths. Fig. 7 (c) and (d) comJ
=
2
!
1
1.78
J
=
2
!
1
2.20
4.92
(0.10)
3.18(0.25)
43.04
and
[C
I
]
emissionarefoundtowards
thewest.
These structures
spectrallyasmeasuredbyDamianietal.
(2017).
Assuming
optically
thin
(
13,
8)
pareJ theII] velocity integrated intensity map with WISE 3.4
where L is the line luminosity in L , S V is the velocity inte= 4 [C!extend3
1.25= 4 ! 3 for higher
0.74velocities(0.04)veil in2.89(0.10) range
110.50further
towards theJ north-east
in thecomponethis warm
the velocity
of 2[C II] in9.71
grated flux in Jy km s 1 , ⌫ is the transition frequency in GHz and emission from1
6 J = 6 !
5 of 12 km
1.12s 1 to 15 km
J =s6 1!
5 is quite0.28
5.12(0.23) of3.31(0.57)
24.97
range. Thissimilar to ular
the case
–
7
km
s
and
a
kinetic
temperature
500
K(Tielens&
Holcloud
http://archive.nrao.edu/nvas/
DL is the distance in kpc. A total CO luminosity of LCO = 9.5 L(0, 0)
J = 8 !of7M17
1.04
8!
7 map shows
0.10a strongthe
9.94 (0.06)(0.14)
1985b),
theand[C 43
II4.11] column
density 114.49
using 7 NRAO/VLA
13
SW,Survey,
where(c)the2005-2007[C JII=] channel
spatial[C II]
Archive
AUI/NRAO
Fig. 4. Colour maps of the integrated
intensityfor(leftright) of theJ = 1 !and0 transitionof 12 CO,lenbach
CO andaverageweofcalculated
H↵40,
41,42
lines
was calculated
thetoobserved
transitions
by accounting
17
2
Eq.Dec.(J2000)
A.1 with a =peak24value
at an o↵setofm00 ⇠ 9.6 ⇥ 10
3.88
(0.16)
4.01cm(0.36)
37.66
association with [C I] and CO channel maps only at intermediate
toward Her 36 which is the central( ↵ =of0,the missing= 0) at transitions.
R.A.(J2000)A=[C18Ih] 0340.3
220 43of, marked withanfor theposition
luminosities609µms and00
2)
a Calculated for di↵erent
1 of 18 gradients↵ = 30(30,,= 1000emission)9.72
from(0.07)
Her
36.inTaking
C/H⇠ 1.2 ⇥100.10
10 4 , Article
with theirWhileerrorsatinlower
the plot
number,page
4.86 1.
(0.16)10 km
s 10
to 24 km s 1 velocities.
(<10askm s 1 )asterisk. The contour levels are10%(> 3 ⇥ rms,
given
100%
in steps whichof 10%isofa the( corresponding
peak
given
Table
obtained,line
luminosityof L
0.11 1)
L towas[C I] in=Table
assuming(60,all carbon
ionized form
and(0.09)
from the column
in Fig.
11. (>24 km s 1 ) velocity channels, [C24
All maps are plotted using original
bottom leftsinceof each27) being
10.22in(0.04)
3.34
112.83 indicated andhigherII] emission islowerbeam
limit sizes
to theshown
total in[Cthe
I] luminosity
themap.
[C I] 370 µm
202
density of H2 , N(H2 ) = 9.4 ⇥ 10 cm (Av /mag) (Kau↵mann
mostly not associated with the other tracers of dense and di↵use
line was not observed. The estimated [C II] luminosity is L[C II]
et al. a2008),
we derivedposition
the visualextinction. A maximum
The referenceis thatof HerAv36.
Table 1. Line parameters of observedgas (Pérez-Beaupuits et al. 2015c). Notably, our [C II] channel= 95.8 Ltransitions.
. Similar to the mass estimations in Sect. 4.1, these extinctionof Av ⇠ 4.25 is obtained at the same position whereb
Table 2. 12 CO, 13 CO, [C I] and [C II] line parameters.

12

h

5

m

s

0

00

018 March 22; in original form 2017 October 13

304 Tes$ng the Larson rela$ons in massive clumps
Traﬁcante,1‹ A. Duarte-Cabral,2 D. Elia,1 G. A. Fuller,3 M. Merello,1 S. Molinari,1 N. Pereqo,2 E. Schisano1 and A. Di Giorgio1#56 A.
MNRAS 477, 2220

ABSTRACT

We tested the validity of the three Larson relations in a sample of 213 massive clumps
selected from the Herschel infrared Galactic Plane (Hi-GAL) survey, also using data from
the Millimetre Astronomy Legacy Team 90 GHz (MALT90) survey of 3-mm emission lines.
The clumps are divided into five evolutionary stages so that we can also discuss the Larson
relations as a function of evolution. We show that this ensemble does not follow the three Larson
relations, regardless of the clump’s evolutionary phase. A consequence of this breakdown isthat the dependence of the virial parameter α vir on mass (and radius) is only a function of
the gravitational energy, independent of the kinetic energy of the system; thus, α vir is not agood descriptor of clump dynamics. Our results suggest that clumps with clear signatures of
infall motions are statistically indistinguishable from clumps with no such signatures. The
observed non-thermal motions are not necessarily ascribed to turbulence acting to sustain thegravity, but they might be a result of the gravitational collapse at the clump scales. This seems
to be particularly true for the most massive (M ≥ 1000 M⊙ ) clumps in the sample, where
exceptionally high magnetic fields might not be enough to stabilize the collapse.
Key words: surveys – stars: formation – stars: kinematics and dynamics – stars: massive –stars: statistics – infrared: stars.

Larson (1981)
0.5(ii) power-law
Clouds are in approximately
with a virial
1. A size–linewidth
relaQon, σvirial∝ Requilibrium,
2/Evirial
R/GM ≃ 1, where
M isathemass parameterofparameter α vir = Ekin
G = 5σ equilibrium,
2. Clouds are approximately
withvirialdominated by highly supersonic
the region
and G is the gravitational constant. This relation implies
2R/GM
=
E
/E
=
5σ
1 of the system Ek ∝ σ 2 M is of similar intensity
αvircloudsk G that the kinetic≃
persions in giant molecularenergyas the gravitationalenergy ofntheon et al. 1987), massive
clumps density
3. A volume
n-size relaQon,∝system,
R−1.1 EG ∝ M2 /R.

. Urquhart et al. 2014; Traficantemassive cores (size ≃ 0.1 pc; Zin-(iii) A volume density n-size relation, n ∝ R−1.1 . This relation
implies that GMCs are universal structures, with a mostly uniform

25

Table 2. Classification scheme of our 213 clumps following the evolutionary scenario described in Merello et al. (in preparation). In column 2, we give thewavelengths at which each evolutionary stage becomes bright. In column 3, we list the surveys that have a visible counterpart in the data. Column 4 shows themedian value of L/M and column 5 gives the number of identified objects.
Evolutionary
phase
Starless
Protostellar MIR dark
Protostellar MIR bright
YSOs
H II regions (UCH II + ext. H II)

Brightwavelengths (µm)

Survey>70
Hi-GAL≥70
Hi-GAL≥21
Hi-GAL, MIPSGAL, WISE
≥8
Hi-GAL, MIPSGAL, WISE, RMS
relationsmassiveclumps≥8 and radio emission
Hi-GAL,WISE,
RMS 2227
Larsonrelationsin massive
clumps2225
2228 Larson
A.
Traficante
etinMIPSGAL,
al.

L/M(L⊙ /M⊙ )
1.7
8.4
6.4
32.1
36.3

Count
14
12
106
25
56 (42+14)

sources. In other words, the concept of virial equilibrium can bemisleading: in a sample that violates the three Larson relations, thevirial parameter varies independently of the kinetic energy of thesources.
Note that these results are valid for an ensemble of clumps, where
the non-thermal motions are estimated with a single gas tracer.
However, this common choice biases the observations towards regions with similar volume densities within each clump, regardlessof the physical properties of the clump (i.e. mass and size). The
correlation between α vir and mass seems to disappear when different surveys of clumps and cores observed with different tracers
are combined (Kauffmann et al. 2013). These results might differ
from the analysis of the energy balance within single star-forming
regions. They might all be near virial equilibrium, as predicted by,for example, Lee & Hennebelle (2016), but the kinetic energy inthe more massive clumps might not be properly measured.
The virial parameter determined for an ensemble of massive
clumps might not be a good descriptor of clump dynamics, which
Figure 4. First Larson relation: velocity dispersion σ as a function of radius is a hypothesis that we explore in the next section.
Figure 5. The first Larson relation distribution for our clumps, divided for
R. The dark grey dashed line is the original Larson relation, σ ∝ R0.38 , anddifferent
evolutionary
The relation
states thatquantity
σ /R0.5
9. Virial
parameterstages.
distribution
as a functionof themass.
The redthe light grey dash-dotted line is the revised Heyer & Brunt (2004) relation, 6FigureIN
FA LbeL Ia NG P ROP E RT
Idistributions
ESshouldconstant.
Instead,
the
span
a
range
of
values0.56
Figureσ2226
8.
Mass–radius
distribution.
The
line
black
dashedlineisthe∝ R . Theis weak,etwith
A.correlation
Traficanteal.a Pearson coefficient of ρ = 0.26. dashed line is the best fit to our data, which gives a slope α = −0.56 ± 0.04. foreach
evolutionary
If thevirialparameterphase.
is not a good descriptor of the energy balance
Kauffmann Figure& Pillai(2010)
star formation
2. L/M
ratiomassive
distributionof our 213threshold,clumps divided
intothevarious1.33
inan
ensemblebeM(r) >
M⊙first(R/pc),phases.
andthe greendash-dottedline is spans
the revised
evolutionaryTheentire sample
more than four orders ofof massive clumps, then these regions−1might
4.1870The
Larsonrelation:
linewidth–size−0.5
median values[2.32,2.21,2.66,2.73,2.60]
kmofsthepc
in theboundofandcollapsing independently
valueof
1.42 . The red
Baldeschi etmagnitude
al. (2017) threshold,M(r)is >an1282
M⊙ (R/pc)
and there
evident
increaseof the L/M gravitationally
ratio going from
One wayprotostar
to exploreMIR
this dark,
hypothesis
is to MIR
comparethe YSOs
properties
protostarbright,and H II
The first Larson relation shows the proportionality
between the size α vir . starless,
2.38
±
0.10
dashed line starlessis the besttofitH
to IIourregions.
sample, M ∝ R.
show signsphases,of infallmotions withclumps orthathave forregionsthat
evolutionary
respectively.
Altogether
dividedof GMCs and the non-thermal motions of the gas in the region of clumps
signatures.
different
evolutionary phases, these results suggest that the first(Larson 1981). This relation has often been considered to be the no such+
3. theGalacticof breakthe 213 clumps investigated in this work.
(1–0) transition,
The
MALT90surveyalsoFigure
observed
Larsonrelationtypically
breaksdown
atHCO
clumpdistribution
scales,and this
Protostar
MIR
which
when
a clump becomesvisibleof(ii)
interstellar
turbulence.
The
interstellar
mediummodelled
Theresultmass–radius
diagramis alsodark,a useful
tool isto
investigate
an opticallythick line
anddifferent
aColours
good internal
tracerofconditions
infall motions(Fuller,and symbols
represent
clumpsat different evolutionary stages. Overis
notbecause
of
the
of
these
objects.
asaturbulent
fluiddominatedbyshocksfollowsapower-spectrumat are70 µm
is stilldark instars.
the TheMIR,
the emission
is too faint toclumps thatlikelybut
to itformhigh-mass
vastormajorWilliams & Sridharan 2005;
al. 2013;
He et
al.
2015;
laidRygl
isrelation
theetGalactic
model
Hou
al. (2009),
the same model used to discussThe observed first Larson
impliesthat of
one,
or et
both,
of
a Burgers-likepower
spectrum
relationofthe form RThese∝ σ 0.5 (i.e.
beidentified.
clumpsarebrightat
all
Hi-GALwavelengths,
et al. 2017). Thesemotionscanbe identified
by looking
ity ofE theseclumpscan form high-mass2007),
objects,
following
therepro- Traficante
the
clump
distribution
inthisElia
et al. (2017).
All clumps are in the IV Quadthe
other
two
relations
must
not
be
followedbyensembleofk−2 ; noe.g. counterparts
McKee & Ostrikerascenariothatk ∝withat bluesurveys.
asymmetries in the spectra corresponding to single-peakedin the MIPSGAL,
WISE& and RMS
Figure 11.
panel:
α virandversusmass diagram for the 21 clumps that
empirical
thresholds determinedby Kauffmannrant, mostly distributed across the Crux-Scutum
andUpperNormaarms
in theclumps.
ducesmass–radiuslarge-scale observations(Padoan & Nordlund
2002; McKee(1–0)
avoids the risk that asymmetries in the
N2 H+
(iii)
Protostar
MIR
bright,
which
is
when
the
clumpsalsospectra,
be- which inter-arm
showkpc
evidenceof infalling
motions inoftheir HCO+ (1–0) spectra. Lower
region. The red-filled dot at [0,8.5]
represents
the position
Pillai&(2010)
and2007).
Baldeschi
et al.the(2017).
but totwobreakclumpsOstriker
However,
relationAll
seems
in mas- HCO+ (1–0) spectra might be
the result of obscuration by surroundpanel: same diagram for the remaining 192 clumps of the sample. The
comevisibleMIR clouds.
andataredifferent
theirbolometricluminosity
increases
Figure
6. Histogram
of in
the
α
evolutionary
stages.
vir distribution
thedistribution
Sun. The ofblack
circle separated
at [0,0] kpc
is theslopes
Galactic
Centre.
Figure 7.(Chira
Surfacedensity
the clumps
for eachsive
clumpsembeddedin themolecularFor
example,Caselli&(G333.449–00.183
and G338.917+00.382)
above
the former,
etLarsonal.
2014).
of the best linear fit in the two plots are the same as the uncertainties,
The majority of the clumps have α vir ≤ 2, independent of the evolutionary ing filaments
4.2
The
secondrelation:
the virial
equilibrium
evolutionary
phase.
The blackvertical lines
correspondto the median values
significantly.
These
clumps
have
at
least
one
counterpart
in
one
of
phase. The
blackdashed line
toB
α vir
= 1 and star-forming
the grey dash(1995)
observed
the corresponds
Orion A andG338.917+00.382
high-mass
Weofidentified
eye all clumps
withinfallsignatures" in eachby
evolutionary
stage.
Thewell-definedsurface density
valuesspan more than suggesting that the clump dynamics is independent of the values of α vir .
and allMyers
but
three
clumps(G333.449–00.183,
anddotted line corresponds to α vir = 2.
one pass
order ofmagnitude, and there is no clear evidence of a trend
theandMIPSGAL,WISE andbetween
MSXsizesurveysbut theynottheis, second
all clumps
withrelation
double-peaked
blue-asymmetric
HCO+ with
regions
they found a correlationand velocity
dis- –dothat
The
Larsonstates that
GMCs are approximately
evolution.
G343.938+00.097)
are above the0.21latter, which is a more stringent(1–0)
spectra
–
and
we
found
21
clumpswith
suchproperties.
, significantly
lower
than
that foundet al. 2013).
been in- of the virial parameter distributions in the two subsamples are thepersionof theform 3.
R to∝
σ classified
in virial equilibrium. The virial parameter α vir has often Thisin our
Table
Range
of
values ofasα vir
RMS
criteria
be
YSOs(Lumsdenclumps for each stage of evolution.
threshold.
massiveclumpshave
shown
thatequilibrium
spanthatalmostordersnumberlikely as
represents
a subsample
oftheyclumpshavetwoparsecand EGof same.
in GMCs.
Similarresults
have beenobtainedthe survey ofhave
high-reachedterpreted
representativeofthebetween
E
k
(iv)
YSOs,
which
is whentheinprotostars
thezeromagnitude in surface densities (Urquhart et al. 2014; Traficante+mass diagrams for the 21 clumps
Fig. 11 showstheisα virnoversus
scalewheninfall2015a;
motions.
Manysuchclumpspresent
red-asymmetricTomassconcludethisEvolutionary
section,we haveshownetthat
therelationsthefindings
MerelloetHCO
al. (inthere
signifiphase
α virLarson
otherforcesasmagnetic
fieldsarenotinvolved(andpreparation),star-forming
regionsof
Shirleyal.
(2003).
At
the
same
time,
et al. all
Svobodaet al.
2016;
Eliaofet al.
2017).
+
age main sequence and have also become bright in the
NIR
regime.
spectrawithblue-asymmetric
HCO
(1–0)
spectradespitethepresenceof
infallmotions
that
can
be
The
surfacedensities
of
the
213
clumpsanalysed
in
this
work
are
Starless
0.07
≤
α
≤
7.98
virfound nodifference
ininL/M
betweenMIR dark and MIR bright sources, (upper panel) and for the redo nototherdescribe
theMIR
properties 0.13
ofobjectsan≤ ensembleof correlation(masassuming
spherical
andcant
homogeneousdensity
distribution;
Bertoldi
surveysofdynamicalmassive star-forming
shown
in Fig.
7. The surfacespectra
densities
the range tracers
0.13 ≤ " ≤ maining 192 clumps (lower panel), with the respective fits. The
dark
αthe
1.71
They areclassified
as YSOs insurvey. identified
vir ≤RMS
withblue-asymmetricin arehigher-density−2,
spanning
more
than
one
order
of
magnitude
in
each
8.57
g
cm
=
α
=
1
or,
&
McKee
1992).
The
virialequilibrium
impliesα(Plumeetal.
1997;Ballesteros-Paredes
et
al.
2011;
Traficante
et
al.
MIR
bright
0.07
≤
α
≤
3.52
vir
suggestingthat
thepresenceofa
MIR source in a clump does notvireqsive) clumps. (v)In theYSOs
nextregions,section, we
exploreimplications of(Wyrowski
et al.of2016).
infallingclumpsmight
not always
0.20≤thermal
α vir ≤ 2.03bremsstrahlung
where
thesomeemission
evolutionary
phase. Also,
We foundmedianvaluesof [0.83,
0.45, 0.80, linear fits in the diagrams in Fig. 11 have very similar slopes within
ifa collapsing
cloud
modelledprotostar
as an alter
isothermal(Bonnor–Ebert)2018) or evenH
anIIIIregions
inverse correlation(Wu
al. ≤2010)−2 is significantlyH
0.05≤etα vir
12.81between size
the
total
luminosity
ofuncertainties,
the cold dustenvelope.
instarless,
24-µm
dark,
protostar
240.98,
0.99]
g
cm
the
suggestingthat the observed values of the virial
this evidence.
showstars
the expected
the gas ionized in the envelope of the more massive
can beblue-asymmetric line profiles (Smith et al. 2013).

26

Dipartimento di Fisica e Chimica, Università di Palermo, Via Archirafi 36, 90123 Palermo, Italy

304 Structure of X-ray emiJng jets close to the launching site: from embedded to disk#57 bearing sources

Received 1 December 2017; accepted 15 March 2018

S. Ustamujic1, S. Orlando2, R. Bonito2, 3, M. Miceli3, 2, and A. I. Gómez de Castro1
ABSTRACT
A&Ap in pressContext. Several observations of stellar jets show evidence of X-ray emitting shocks close to the launching site. In some cases, theshocked features appear to be stationary, also for Young Stellar Objects (YSOs) at di↵erent stages of evolution. We study the case of
HH 154, the jet originating from the embedded binary Class 0/I protostar IRS 5, and the case of the jet associated to DG Tau, a more
evolved Class II disk-bearing source or Classical T Tauri star (CTTS), both located in the Taurus star-forming region.
Aims. We aim at investigating the e↵ect of perturbations in X-ray emitting stationary shocks in stellar jets; the stability and detectability in X-rays of these shocks; and explore the di↵erences in jets from Class 0 to Class II sources.
Methods. We performed a set of 2.5-dimensional magnetohydrodynamic numerical simulations that modelled supersonic jets ramming into a magnetized medium. The jet is formed by two components: a continously driven component that forms a quasi-stationary
shock at the base of the jet; and a pulsed component constituted by blobs perturbing the shock. We explored di↵erent parameters forboth components. We studied two cases: a jet less dense than the ambient medium (light jet), representing the case of HH 154; anda jet denser than the ambient (heavy jet), associated with DG Tau. We synthesized the count rate from the simulations and comparedwith available Chandra observations.
Results. Our model explains the formation of X-ray emitting quasi-stationary shocks observed at the base of jets in a natural way,
being able to reproduce the observed jet properties at di↵erent evolutionary phases (in particular, for HH 154 and DG Tau). The jet iscollimated by the magnetic field forming a quasi-stationary shock at the base which emits in X-rays even when perturbations formed
by a train of blobs are present. We found similar collimation mechanisms dominating in both heavy and light jets.
Conclusions. We derived the pysical parameters that can give rise to X-ray emission consistent with observations of HH 154 and
DG Tau. We have also performed a wide exploration of the parameter space characterizing the model; this can represent a useful tool
to study and diagnose the physical properties of YSO jets over a broad range of physical conditions, from embedded to disk-bearing
sources. We have shown that luminosity did not change significantly in variable jet models for the range of parameters explored.
Finally, we provided an estimation of the maximum perturbations that can be present in HH 154 and DG Tau taking into account theavailable X-ray observations.
Key words. ISM: jets and outflows – magnetohydrodynamics (MHD) – X-rays: ISM – stars: pre-main sequence – ISM: individual
objects: HH 154 – stars: individual: DG Tau

Introductionung Stellar Objects (YSOs) are stars at their early stages ofdisk, and most of its complex phenomenology can be mode
as a star interacting with a circumstellar accretion disk
27 (Stone
Norman 1994; Stone et al. 1996; Romanova et al. 2011; Orlan

304 Star forma$on toward the H II region IRAS 10427-6032
Vaidya⁎,a,
A B S Kaushar
T R A
C T Souradeep Bhaqacharyaa,b, Vatsal Panwara,c, Manash R. Samald, Wen-Ping Chend, Devendra K. Ojhae#58
New Astronomy 63, 27

The formation and properties of star clusters formed at the edges of the H II regions are poorly known. In thispaper, we study stellar content, physical conditions, and star formation processes around a relatively unknown
young H II region IRAS 10427-6032, located in the southern outskirts of the Carina Nebula. We study this region
by making use of the near infrared (near-IR) data from VISTA, mid infrared (mid-IR) from Spitzer and WISE, farinfrared (far-IR) from Herschel, sub-mm from ATLASGAL, and 843 MHz radio-continuum data. Using multi-band
photometry, we identify a total of 5 Class I and 29 Class II young stellar object (YSO) candidates, most of which
newly identiﬁed, in the 5′ × 5′ region centered on the IRAS source position. Modelling of the spectral energy
distribution for selected YSO candidates using the radiative transfer models shows that most of these candidates
are intermediate mass YSOs in their early evolutionary stages. A majority of the YSO candidates are found to becoincident with the cold dense clump at the western rim of the H II region. Lyman continuum luminosity calculation using radio emission indicates the spectral type of the ionizing source to be earlier than B0.5–B1. We
identiﬁed a candidate massive star possibly responsible for the H II region with an estimated spectral type
B0–B0.5. The temperature and column density maps of the region constructed by performing pixel-wise modiﬁed
blackbody ﬁts to the thermal dust emission using the far-IR data from the Herschel show a high column density
shell-like morphology around the H II region, and low column density (0.6 × 1022 cm −2 ) and high temperature(∼ 21 K) matter within the H II region. Based on the morphology of the region in the ionized and the moleculargas, and the comparison between the estimated timescales of the H II region and the YSO candidates in theclump, we argue that the enhanced star-formation at the western rim of the H II region is likely due to compression by the ionized gas.

dense clumps in which star formation is likely to occur (Kendrew et al.,
28 the2016). An expanding H II region may trigger star formation via