Herschel survey and modelling of externally-illuminated photoevaporating protoplanetary disks

TL;DR

This study combines Herschel and ALMA observations with a new PDR model to analyze externally-illuminated photoevaporating disks, revealing their energetic regimes, physical parameters, and self-regulated mass-loss rates, advancing understanding of disk evolution.

Contribution

The paper introduces a 1D PDR model tailored for proplyds, linking observations with physical parameters and identifying a unique energetic regime affecting disk photoevaporation.

Findings

Derived densities at disk surfaces (~10^6 cm^-3)

Estimated mass-loss rates (~10^-7 M_sun/yr)

Identified a self-regulated energetic regime

Abstract

Protoplanetary disks undergo substantial mass-loss by photoevaporation, a mechanism which is crucial to their dynamical evolution. However, the processes regulating the gas energetics have not been well constrained by observations so far. We aim at studying the processes involved in disk photoevaporation when it is driven by far-UV photons. We present a unique Herschel survey and new ALMA observations of four externally-illuminated photoevaporating disks (a.k.a. proplyds). For the analysis of these data, we developed a 1D model of the photodissociation region (PDR) of a proplyd, based on the Meudon PDR code and computed the far infrared line emission. We successfully reproduce most of the observations and derive key physical parameters, i.e. densities at the disk surface of about $10^{6}$ cm$^{-3}$ and local gas temperatures of about 1000 K. Our modelling suggests that all studied disks are found in a transitional regime resulting from the interplay between several heating and cooling processes that we identify. These differ from those dominating in classical PDRs, i.e. grain photo-electric effect and cooling by [OI] and [CII] FIR lines. This energetic regime is associated to an equilibrium dynamical point of the photoevaporation flow: the mass-loss rate is self-regulated to set the envelope column density at a value that maintains the temperature at the disk surface around 1000 K. From our best-fit models, we estimate mass-loss rates - of the order of $10^{-7}$ $\mathrm{M}_\odot$/yr - that are in agreement with earlier spectroscopic observation of ionised gas tracers. This holds only if we assume an evaporation flow launched from the disk surface at sound speed (supercritical regime). We have identified the energetic regime regulating FUV-photoevaporation in proplyds. This regime could be implemented into models of the dynamical evolution of protoplanetary disks.

Tables

Table 1. Table 1: General properties of the studied proplyds.

Proplyd name	RA	DEC	Apparent	Apparent	FUV field $G_0$	Heliocentric
(parent nebula)	(J2000.0)	(J2000.0)	envelope size	disk size	(Habing field)	distance (parsec)
105-600 (Carina)	10 46 32.97a𝑎aa𝑎a$a$Sahai et al. (2012)	$-$ 60 03 53.50a𝑎aa𝑎a$a$Sahai et al. (2012)	$9.5″\times 3.7″$b𝑏bb𝑏b$b$Vicente (2009b)	…	$2.2\times {10}^4$b𝑏bb𝑏b$b$Vicente (2009b)	$2300$c𝑐cc𝑐c$c$Smith et al. (2003)
HST10   (Orion)	05 35 18.22d𝑑dd𝑑d$d$Ricci et al. (2008)	$-$ 05 24 13.45d𝑑dd𝑑d$d$Ricci et al. (2008)	$2.6″\times 1.0″$e𝑒ee𝑒e$e$Chen et al. (1998)	$0.4″\times 0.1″$e𝑒ee𝑒e$e$Chen et al. (1998)	$2.4\times {10}^5$f𝑓ff𝑓f$f$Störzer & Hollenbach (1998)	$414$g𝑔gg𝑔g$g$Menten et al. (2007)
244-440 (Orion)	05 35 24.38d𝑑dd𝑑d$d$Ricci et al. (2008)	$-$ 05 24 39.74d𝑑dd𝑑d$d$Ricci et al. (2008)	$5.6″$hℎhhℎh$h$Vicente & Alves (2005)	$0.86″\times 0.69″$hℎhhℎh$h$Vicente & Alves (2005)	$1\times {10}^5$b𝑏bb𝑏b$b$Vicente (2009b)	$414$g𝑔gg𝑔g$g$Menten et al. (2007)
203-506 (Orion)	05 35 20.32d𝑑dd𝑑d$d$Ricci et al. (2008)	$-$ 05 25 05.55d𝑑dd𝑑d$d$Ricci et al. (2008)	…	$0.75″\times 0.61″$j𝑗jj𝑗j$j$Noel (2003)	$2\times {10}^4$i𝑖ii𝑖i$i$mean value of Marconi et al. (1998) and Walmsley et al. (2000)	$414$g𝑔gg𝑔g$g$Menten et al. (2007)

Table 2. Table 2: Observing strategy for Herschel spectroscopy.

Line	Target	Instrument	Observing mode	Band	Beam size
					(FWHM or spaxel witdth)
[CII] 158 µm	105-600	HIFI - HRS	Dual beam switch	7b (1795.2 - 1902.8 GHz or 158 - 177 µm)	11.3 ″
	HST10		Load chop
	244-440		Load chop
CO (7-6)	105-600	HIFI - HRS	Dual beam switch	3a (807.1 - 851.9 GHz or 352 - 372 µm)	28.4 ″
	HST10		Frequency switch
	244-440		Frequency switch
CO (10-9)	105-600	HIFI - HRS	Frequency switch	5a (1116.2 - 1240.8 GHz or 242 - 269 µm)	18.6 ″
	HST10
	244-440
CO (15-14), (17-16), (19-18)	105-600	PACS	Chopping nodding	51-73 µm and 103-220 µm	9.4 ″
	HST10		Unchopped grating scan
	244-440		Unchopped grating scan
[OI] 63 µm	105-600	PACS	Chopping nodding	51-73 µm and 103-220 µm	9.4 ″
	HST10		Unchopped grating scan
	244-440		Unchopped grating scan

Table 3. Table 3: Features of the observed lines ( Herschel plus published/archival data available).

Target	Line	Wavelength	Instrument	Integrated intensity			FWHM	$v_{\mathrm{LSR}}$
		(µm)		(W m-2)		(W m-2 sr-1)	(km s-1)	(km s-1)
105-600	[CII] 158 µm	157.74	HIFI	2.39e-16	(12%)	3.68e-07c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	6.4	-23.1
	CO (7-6)	371.65	HIFI	1.18e-17	(11%)	…	2.7	-22.7
	CO (10-9)	260.24	HIFI	2.26e-17	(14%)	…	2.5	-22.8
	CO (15-14)	173.63	PACS	5.50e-18	(24%)	…	…	…
	CO (17-16)	153.27	PACS	7.59e-18	(24%)	…	…	…
	CO (19-18)	137.20	PACS	¡ 2.65e-17a𝑎aa𝑎a$a$Upper limit determined by the noise (peak of the line supposed to be smaller than three times the noise level).	(20%)	…	…	…
	[OI] 63 µm	63.18	PACS	2.90e-15	(12%)	4.47e-06c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	…	…
	CO (3-2)	866.96	APEX/FLASH+	2.98e-18	(20%)	…	3.2	-22.7
	CO (4-3)	650.25	APEX/FLASH+	6.67e-18	(20%)	…	3.1	-22.7
	CO (6-5)	433.56	APEX/CHAMP+	1.26e-17	(20%)	…	3.5	-22.6
	CO (7-6)	371.65	APEX/CHAMP+	1.40e-17	(20%)	…	4.1	-22.3
	13CO (3-2)	906.83	APEX/FLASH+	7.12e-19	(20%)	…	2.2	-22.6
	HCO+ (4-3)	840.41	APEX/FLASH+	2.97e-19	(20%)	…	2.0	-22.6
	HCN (4-3)	845.66	APEX/FLASH+	1.72e-19	(20%)	…	1.8	-22.8
HST10	[CII] 158 µm	157.74	HIFI	¡ 2.27e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(12%)	4.72e-06b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	3.9	21.3
	CO (7-6)	371.65	HIFI	¡ 4.32e-18a𝑎aa𝑎a$a$Upper limit determined by the noise (peak of the line supposed to be smaller than three times the noise level).	(11%)	5.84e-06a,c𝑎𝑐a,ca,c𝑎𝑐a,cfootnotemark: $a,c$	…	…
	CO (10-9)	260.24	HIFI	¡ 2.63e-17a𝑎aa𝑎a$a$Upper limit determined by the noise (peak of the line supposed to be smaller than three times the noise level).	(14%)	3.56e-05a,c𝑎𝑐a,ca,c𝑎𝑐a,cfootnotemark: $a,c$	…	…
	CO (15-14)	173.63	PACS	¡ 7.15e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	9.68e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	CO (17-16)	153.27	PACS	¡ 3.47e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	4.70e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	CO (19-18)	137.20	PACS	¡ 1.44e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	1.95e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	[OI] 63 µm	63.18	PACS	¡ 1.01e-13b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(12%)	2.10e-03b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	H2 1-0 S(1)	2.12	HST/NICMOS	1.85e-18c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	(20%)	2.50e-06	…	…
	[OI] 6300Å	0.63	HST/WFPC2	1.11e-17c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	(34%)	1.50e-05	…	…
244-440	[CII] 158 µm	157.74	HIFI	¡ 4.68e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(12%)	8.09e-07b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	2.2	21.7
	CO (7-6)	371.65	HIFI	¡ 2.69e-18a𝑎aa𝑎a$a$Upper limit determined by the noise (peak of the line supposed to be smaller than three times the noise level).	(11%)	2.46e-07a,c𝑎𝑐a,ca,c𝑎𝑐a,cfootnotemark: $a,c$	…	…
	CO (10-9)	260.24	HIFI	¡ 1.53e-17a𝑎aa𝑎a$a$Upper limit determined by the noise (peak of the line supposed to be smaller than three times the noise level).	(14%)	1.40e-06a,c𝑎𝑐a,ca,c𝑎𝑐a,cfootnotemark: $a,c$	…	…
	CO (15-14)	173.63	PACS	¡ 4.19e-15b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	3.83e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	CO (17-16)	153.27	PACS	¡ 2.72e-15b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	2.49e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	CO (19-18)	137.20	PACS	¡ 8.16e-16b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(13%)	7.45e-05b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	[OI] 63 µm	63.18	PACS	¡ 2.06e-13b𝑏bb𝑏b$b$Upper limit determined by the observed background emission.	(12%)	3.55e-04b,c𝑏𝑐b,cb,c𝑏𝑐b,cfootnotemark: $b,c$	…	…
	OH 84 µm	84.45	PACS	6.83e-18	(13%)	6.24e-07c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	…	…
203-506	HCO+ (4-3)	840.41	ALMA	3.04e-20c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	(10%)	3.60e-09	4.0	9.5
	H2 1-0 S(1)	2.12	CFHT	7.01e-18c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	(11%)	8.30e-07	…	…
	[OI] 6300Å	0.63	HST/WFPC2 - MUSE	1.01e-16c𝑐cc𝑐c$c$Values converted using the apparent size given in Table 1 and assuming that all the lines are coming from the disk except for [CII] 158 µm and [OI] 63 µm lines that are supposed to originate in the envelope.	…	1.20e-05	…	…

Table 4. Table 4: Continuum emission of the Carina candidate proplyd 105-600.

Wavelength	Instruments	Continuum		AORa𝑎aa𝑎a$a$Astronomical observation request. or
(µm)		emission (mJy)		reference
3.55	IRAC	3.5	(11 %)	12914432
4.49	IRAC	3.5	(11 %)	12914432
5.73	IRAC	13.9	(11 %)	12914432
7.87	IRAC	28.3	(11 %)	12914432
24.00	MIPS	¡ 1144.2b𝑏bb𝑏b$b$Values in MJy/sr.		23788800
57.95	PACS	1153.5	(36 %)	OT2_oberne_4
63.22	PACS	839.4	(45 %)	OT2_oberne_4
68.70	PACS	1013.1	(25 %)	OT2_oberne_4
70.00	PACS	1483.0	(11 %)	1342255062
137.25	PACS	1086.4	(15 %)	OT2_oberne_4
153.19	PACS	1144.2	(16 %)	OT2_oberne_4
160.00	PACS	1035.0	(11 %)	1342255062
173.75	PACS	1122.7	(21 %)	OT2_oberne_4
189.50	PACS	1007.7	(43 %)	OT2_oberne_4
250.00	SPIRE	50.3b𝑏bb𝑏b$b$Values in MJy/sr.	(14 %)	1342255061
350.00	SPIRE	10.9b𝑏bb𝑏b$b$Values in MJy/sr.	(36 %)	1342255061
350.00	LABOCA	196.0	(34 %)	Sahai et al. (2012)
870.00	SABOCA	¡ 40.0		Sahai et al. (2012)

Table 5. Table 5: Parameters of the models for the four studied proplyds.

Symbol	Description	105-600	HST10	244-440	203-506
$a_{\mathrm{env}}$	Envelope semi-major axis	$4.75″$	$1.3″$	$2.8″$	…
$b_{\mathrm{env}}$	Envelope semi-minor axis	$1.85″$	$0.5″$	$2.8″$	…
$n_{\mathrm{env}}$	Envelope density	$8.5\times {10}^3$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$2.5\times {10}^5$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$1\times {10}^5$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$4.5\times {10}^4$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.
$T_{D,\mathrm{env}}$	Mean dust temperature in the envelope	$39.8$ Ka𝑎aa𝑎a$a$Values obtained from the best-fit model.	…	…	…
$f{}_C^{\mathrm{PAH}}$	Fraction of gas phase C locked in PAHs	$0.24$ %a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$0.08$ %b𝑏bb𝑏b$b$Value derived by Vicente et al. (2013).	…	…
$S_{\mathrm{disk}}$	Apparent surface of the disk	$3.82$ arcsec2a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$0.0314$ arcsec2	$0.47$ arcsec2	$0.36$ arcsec2
$n_{\mathrm{disk}}$	Mean disk surface layer density	$7\times {10}^5$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$4\times {10}^6$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$2\times {10}^6$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.	$1\times {10}^6$ cm-3a𝑎aa𝑎a$a$Values obtained from the best-fit model.
	Mean dust temperature in the disk	$19.5$ Ka𝑎aa𝑎a$a$Values obtained from the best-fit model.	…	…	…
$G_0$	Radiation field	$2.2\times {10}^4$	$2.4\times {10}^5$	$1.0\times {10}^5$	$2.0\times {10}^4$

Table 6. Table 6: Default values of some parameters used in the Meudon PDR code.

Parameter	Value
Cosmic ray ionisation rate	5 $\times$ 10-17 s-1
Turbulent velocity (Doppler broadening)	2 km s-1
R${}_V{}^{}=A_V/E(B-V)$	3.10
N${}_H{}^{}/E(B-V)$	5.8 $\times$ 1021 cm-2
Dust-to-gas mass ratio $\delta$	0.01
Power-law index of the grain size distribution	3.50
Minimum grain radius $a_{\min }$	3 nm
Maximum grain radius $a_{\max }$	300 nm

Table 7. Table 7: Parameters derived from the fit of the spectral energy distribution for 105-600.

Parameters	Best value	Interval at 1\textsigma
$M_{\mathrm{disk}}$ (M☉)	0.67	0.46 - 0.97
$T_{D,\mathrm{disk}}$ (K)	19.5	17.9 - 21.4
$M_{\mathrm{env}}$ (M☉)	0.048	0.036 - 0.062
$T_{D,\mathrm{env}}$ (K)	39.8	37.7 - 42.4
$f{}_C^{\mathrm{PAH}}$ (%)	0.24	0.22 - 0.25

Table 8. Table 8: Physical parameters estimated with the model for the studied prolyds.

Object	Ptrans	Gas temperature at	Nenv	rdisk	$\dot{M}$ supercritical	$\dot{M}$ subcritical	$\dot{M}$ observed
		the disk surface (K)	(cm-2)	(AU)	(M⊙/yr)	(M⊙/yr)	(M⊙/yr)
105-600	1.00	1100	1.4e+21	2500a𝑎aa𝑎a$a$Disk radius roughly approximated by this study (see Sect. 4.1.2).	7.5e-05	…	$1$e-06
HST10	1.05	1150	2.0e+21	85b𝑏bb𝑏b$b$Disk radius estimated by spatially resolved observation: HST10 from Chen et al. (1998), 244-440 and 203-506 from Vicente & Alves (2005)	5.0e-07	1.7e-08	$\left(7\pm 5\right)$e-07
244-440	1.01	1270	1.9e+21	193b𝑏bb𝑏b$b$Disk radius estimated by spatially resolved observation: HST10 from Chen et al. (1998), 244-440 and 203-506 from Vicente & Alves (2005)	1.4e-06	8.2e-09	$\left(15\pm 7\right)$e-07
203-506	1.12	1080	1.0e+21	150b𝑏bb𝑏b$b$Disk radius estimated by spatially resolved observation: HST10 from Chen et al. (1998), 244-440 and 203-506 from Vicente & Alves (2005)	9.4e-08	9.9e-09	…

Equations

$$L_{\nu }\left(T_{\mathrm{D}}\right)=\tau {}_{\nu }{B}_{\nu }(T_{\mathrm{D}}),$$

$$\tau {}_{\nu }=\tau {}_{\nu {}_0}{\left(\frac{\nu }{\nu {}_0}\right)}^{\beta },$$

$$\sigma =\frac{\tau {}_{250}}{N_{\mathrm{H}}}=2.32\pm 0.3\times 10^{-25}\text{leavevmode}{\mathrm{cm}}^2.$$

$$M=N_{\mathrm{H}}\text{leavevmode}\Omega \text{leavevmode}\mu \text{leavevmode}{m}_{\mathrm{H}},$$

$$L_{\nu }\left(T_{\mathrm{D}}\right)=\sigma \left(\frac{M}{\Omega \text{leavevmode}\mu \text{leavevmode}{m}_{\mathrm{H}}}\right){\left(\frac{\nu }{\nu {}_{250}}\right)}^{\beta }{B}_{\nu }(T_{\mathrm{D}}).$$

$$A_{\mathrm{V}}^{\mathrm{jump}}=\frac{N_{\mathrm{env}}}{\kappa }=\frac{a_{\mathrm{env}}\times n_{\mathrm{env}}}{\kappa },$$

$$r_{\mathrm{g}}=\frac{G{M}_{\ast }}{c_{\mathrm{S}}^2},$$

$$c_{\mathrm{S}}=\sqrt{\frac{\gamma \text{leavevmode}{\mathrm{k}}_{\mathrm{B}}\text{leavevmode}T}{\mu \text{leavevmode}{m}_{\mathrm{H}}}}.$$

$$\dot{M}=\mu m_H{n}_{\mathrm{disk}}{c}_{\mathrm{S}}\times 2\pi (r_{\mathrm{disk}}-r_{\mathrm{g}}{)}^2.$$

$${\dot{M}}=C_{0}\leavevmode\nobreak\ N_{C}\leavevmode\nobreak\ \mu\leavevmode\nobreak\ m_{H}\leavevmode\nobreak\ c_{\mathrm{S}}\leavevmode\nobreak\ r_{g}\Bigl{(}{r_{g}\over r_{d}}\Bigr{)}\exp\bigl{[}-r_{g}/2r_{d}\bigr{]}\,,$$

$$T_{\mathrm{env},\mathrm{IF}}\approx \frac{n_{\mathrm{e}}{T}_{\mathrm{HII}}}{n_{\mathrm{env}}},$$

$$N_{\mathrm{u}}=\frac{8\pi k_{\mathrm{B}}{\nu }^2}{h{c}^3{A}_{\mathrm{ul}}}\int T_{\mathrm{B}}\left({}^{13}\mathrm{CO}\right)d\nu ,$$

$$N_{\mathrm{u}}=\frac{N}{Z}{g}_{\mathrm{u}}\exp \left(-\frac{E_{\mathrm{u}}}{k_{\mathrm{B}}{T}_{\mathrm{ex}}}\right),$$

$$T_{\mathrm{B}}\left({}^{12}\mathrm{CO}\right)=\frac{h\nu }{k_{\mathrm{B}}}\left[\frac{1}{\exp \left(\frac{h\nu }{k_{\mathrm{B}}{T}_{\mathrm{ex}}}\right)-1}-\frac{1}{\exp \left(\frac{h\nu }{k_{\mathrm{B}}{T}_{\mathrm{bb}}}\right)-1}\right],$$

$$M_{\mathrm{mol}}=\mu \text{leavevmode}{m}_{\mathrm{H}}\text{leavevmode}S\text{leavevmode}{N}_{{\mathrm{H}}_2},$$

$$\dot{M}=\mu \text{leavevmode}{m}_{\mathrm{H}}\text{leavevmode}{n}_{\mathrm{e}}\text{leavevmode}S\text{leavevmode}{c}_{\mathrm{S}},$$

Full text

*Aims. *We aim at studying the processes involved in disk photoevaporation when it is driven by far-UV photons (i.e. 6 ¡ E ¡ 13.6 eV).

*Methods. *We present a unique Herschel survey and new ALMA observations of four externally-illuminated photoevaporating disks (a.k.a. proplyds). For the analysis of these data, we developed a 1D model of the photodissociation region (PDR) of a proplyd, based on the Meudon PDR code and we computed the far infrared line emission.

Results. With this model, we successfully reproduce most of the observations and derive key physical parameters, i.e. densities at the disk surface of about $10^{6}$ cm-3 and local gas temperatures of about 1000 K. Our modelling suggests that all studied disks are found in a transitional regime resulting from the interplay between several heating and cooling processes that we identify. These differ from those dominating in classical PDRs i.e. grain photo-electric effect and cooling by [OI] and [CII] FIR lines. This specific energetic regime is associated to an equilibrium dynamical point of the photoevaporation flow: the mass-loss rate is self-regulated to keep the envelope column density at a value that maintains the temperature at the disk surface around 1000 K. From the physical parameters derived from our best-fit models, we estimate mass-loss rates – of the order of $10^{-7}$ M☉/yr – that are in agreement with earlier spectroscopic observation of ionised gas tracers. This holds only if we assume photoevaporation in the supercritical regime where the evaporation flow is launched from the disk surface at sound speed.

*Conclusions. *We have identified the energetic regime regulating FUV-photoevaporation in proplyds. This regime could be implemented into models of the dynamical evolution of protoplanetary disks.

1 Introduction

Protoplanetary disks are found around young stars as a natural consequence of angular momentum conservation in a rotating collapsing core. During a few million years (Zuckerman et al. 1995; Haisch et al. 2001), they are the reservoirs of dust and gas for the formation and the growth of planets. After this relatively short time, only a small fraction of the disk initial mass remains in the form of planets. The main part is supposed to have been accreted onto the central star or physically removed from the disk and one proposed mechanism is photoevaporation (see recent reviews in e.g. Alexander et al. 2014; Gorti et al. 2015). Indeed, disks are subject to highly energetic photons that can sufficiently heat the gas at the disk surface to make it overcomes the gravitational field leading to escape flows and mass-loss.

Photoevaporation of disks caused by Extreme-UV photons (or EUV, $13.6$ eV $<h\nu<0.1$ keV) arising from the central star was first studied by Hollenbach et al. (1994); Yorke & Welz (1996); Richling & Yorke (1997); Clarke et al. (2001) who showed that it can create a gap in the inner part of the disk and erodes it. Gorti & Hollenbach (2009) have investigated the case of Far-UV (or FUV, $6$ eV $<h\nu<13.6$ keV) photons and found that they can cause significant erosion since they can remove the gas from the outer disk where most of the mass is located. Photoevaporation of disks caused by EUV and FUV photons from nearby massive stars has also been analytically and numerically modelled (Johnstone et al. 1998; Richling & Yorke 1998, 2000; Störzer & Hollenbach 1999), in order to understand the morphology and evolution of externally illuminated disks, also known as proplyds (O’dell et al. 1993). These latters generally appear as disks surrounded by a tear-drop-shaped ionised envelope with a bright head facing the main UV source and an elongated tail pointing away. After the work of Johnstone et al. (1998) and Störzer & Hollenbach (1998, 1999), supported by observations (e.g. Bally et al. 1998; Chen et al. 1998), it was suggested that FUV photons penetrate the disk, dissociate molecules and heat the gas to temperatures about 102 to 103 K, thus launching a thermal wind at the disk surface. EUV photons cannot penetrate this flow but ionise the gas on its surface and raise its temperature to 104 K. The ionisation front, delimiting the flow that feeds the envelope of the proplyd, is generally located at several disk radii.

Regions where a strong FUV radiation field drives the physics and chemistry are called PhotoDissociation Regions (PDRs). All of the previously cited models of photoevaporation that include FUV photons are based on PDR models. Tielens & Hollenbach (1985) have developed the first general model describing PDRs for a dense molecular cloud illuminated by a strong FUV radiation field with an intensity $G_{\mathrm{0}}$ between 103 and 106 times the Habing field (Habing 1968), which corresponds to the average FUV flux in the interstellar medium and is valued at $1.6\times 10^{-6}\leavevmode\nobreak\ \mathrm{W\,m^{-2}}$ when integrated between 91.2 and 240 nm. To compute the structure and emission of the cloud, the model includes the main heating mechanisms such as the photoelectric effect, carbon ionisation, H2 dissociation and the main cooling mechanisms by few fine-structure lines of atoms and low-lying rotational lines from molecules for example. Thanks to PDR modelling, Störzer & Hollenbach (1999) have been able to reproduce the observed intensity of H2 ($\nu=1-0$) S(1) and make predictions for the major gas cooling lines in the far-infrared ([OI], [CII], etc.). However so far, these latter lines have never been observed towards proplyds. Walsh et al. (2013) modelled the molecular line emission from a protoplanetary disk illuminated by an external radiation field $G_{\mathrm{0}}=4\times 10^{5}$. Facchini et al. (2016) recently published a study in which they use a 3D PDR model to compute the temperature at the surface of an externally illuminated disk. Their study concerns mainly disks situated in regions of moderate UV fields (up to $G_{\mathrm{0}}=3000$). Those authors also do not compare these models with observations as they are mostly interested in modelling dynamics.

We conducted the first and unique far-IR (FIR) survey of proplyds with the Herschel Space Observatory. The Herschel data were combined with multi-wavelength observations (optical to sub-mm) available in the literature and analysed in the frame of a PDR model specially developed for this case. In this paper the following questions are adressed:

• •

Can the FIR emission of proplyds be modelled as a PDR?

• •

What are the physical properties of this PDR (density, temperature)?

• •

What are the key mechanisms driving the thermal balance (heating and cooling processes) and hence the photoevaporation?

We start by introducing the observations (Sect. 2) and the model (Sect. 3). We then confront them to extract the physical properties of the studied objects (Sect. 4) and use the model to investigate the main mechanisms driving the disk surface properties (Sect. 5) and the consequences on photoevaporation and mass-loss rates (Sect. 6). Sect. 7 contains our concluding remarks.

2 Observations

2.1 Herschel survey

Although it is well established that FUV photons play a key role in the gas heating in PDRs, there is no quantitative study of classical PDR tracers in proplyds so far, mostly due to their small sizes and strong background emission. Since these tracers are mostly found at far-IR wavelengths, and hence have to be observed from space, the ESA Herschel Space Observatory (Pilbratt et al. 2010) was the first to provide a sufficient sensitivity and resolution to observe them.

2.1.1 Targets

To maximise the signal-to-noise ratio for the first survey of this kind, we have carefully selected three proplyds (see Table 1 but note that 203-506 will be presented later) according to their relatively large size (several arcseconds in diameter) and their strong incoming FUV radiation field ($G_{\mathrm{0}}>10^{4}$). These targets are briefly described below:

• •

One of the giant proplyd candidates of Smith et al. (2003) in the Carina nebula (Fig. 1(a)), namely 104632.9-600354 in their paper and hereafter called 105-600 following Mesa-Delgado et al. (2016). It is located amongst the southern pillars at a projected distance of 25.4 arcminutes south-east from \texteta Carinae and at 1 arcminute from the Bochum 11 cluster to the west. It was first studied by Vicente (2009b) and has been observed with HST/ACS in H\textalpha emission (Smith et al. 2010a), revealing a large envelope pointing towards \texteta Carinae and a collimated bipolar jet that confirms the pre-stellar nature of this object. The large size of this object as well as its high molecular content ($M\approx$ 0.35 M☉ , after Sahai et al. 2012) could suggest that the young stellar object is still embedded in the remnant molecular core while the morphology is similar to a proplyd. From this point forward, we will consider that the molecular core inside the object is a disk since the true nature is not important for our study.

• •

HST10 (also known as 182-413, Fig. 1(b)) is the most studied proplyd and was one of the first bright objects in which an embedded disk was seen in silhouette (e.g. O’dell et al. 1993; O’dell & Wong 1996) and in emission through the [OI] 6300 Å and H2 (1-0) S(1) lines (Bally et al. 1998, 2000; Chen et al. 1998). It is located at a projected distance of 56″ from the star \textTheta1 Ori C to the SE, and at 32.6″ from \textTheta2 Ori A to the NW. It is a teardrop shaped proplyd containing a prominent nearly edge-on disk.

• •

244-440 (Fig. 1(c)) is the largest proplyd observed in the Orion nebula cluster. It exhibits a faint nearly edge-on silhouette disk (Vicente & Alves 2005) and a one-sided microjet drived by the central star (Bally et al. 2000). It is located at a projected distance of 142″ from the star \textTheta1 Ori C and at 29″ from \textTheta2 Ori A.

2.1.2 Observing strategy and data reduction

The main concern about observations of proplyds is that they are small objects so their emission is likely to be diluted in the Herschel beam (Fig. 1), and their emission can also be confused with emission from the nebula where they lie. Since the PDR density of proplyds is expected to be high at some locations, especially at the disk surface (about $10^{6}$ cm-3), some PDR tracers such as the [OI] 63 µm and the high-$J$ CO lines ($J=$ 15-14, 17-16 and 19-18) are expected to be bright compared to background nebular emission (which has a lower density) and we thus proposed to observe them with PACS (Poglitsch et al. 2010). With HIFI (de Graauw et al. 2010), the [CII] 158 µm line and fainter lower $J$ CO lines ($J=$ 7-6 and 10-9) were observed assuming that the proplyd emission can be spectrally separated from background/foreground nebular emission with a velocity resolution of $190$ m s-1 or less. Thanks to these two instruments, we were able to observe for the first time in proplyds selected excited CO lines spanning a wide range of densities (critical densities from 105 to 107 cm*-3*, see e.g. Yang et al. 2010), and the [OI] 63 µm and [CII] 158 µm fine-structure lines (critical densities from about 103 to 105 cm*-3*) which are generally the two major cooling lines in PDRs (Tielens & Hollenbach 1985).

From one object or one line to another, different observing modes have been used to avoid contamination by the background emission (see Table 2). For HIFI, the mode ”Dual Beam Switch” uses an OFF position (i.e. away from the source) on the sky if a nearby emission-free zone is expected, as for 105-600 mainly. If there is no such zone, the ”Load Chop” mode allows to use internal loads as references and the mode ”Frequency Switch” make its reference by switching between two frequencies in the local oscillator. In a similar way with PACS, we use the ”Standard Chopping-nodding mode” for 105-600 assuming a vicinity free of emission while the ”Unchopped grating scan” mode was used for the Orion proplyds that enables to go farther away from the nebula to find a clean OFF position.

The data reduction and parts of post-processing were done using the Herschel Interactive Processing Environment (HIPE version 13, see Ott 2010) in order to obtain level-2 calibrated spectra (Fig. 2). Line profiles were then fitted using gaussian functions to extract the line integrated intensities and other characteristics such as the full width at half maximum (FWHM) and the velocity with respect to the local standard-of-rest, $v_{\mathrm{LSR}}$ (see Table 3). Uncertainties have been calculated by quadratically summing the instrumental (see user manuals of the instruments for more details) and fitting uncertainties.

2.1.3 Observing results

Fig. 2 shows the spectra extracted from the Herschel data. A summary of all detected lines and upper limits is given in Table 3.

For the Carina candidate proplyd, 105-600, PACS spectra are not contaminated by the surrounding nebula so that the target emission is directly observed (Fig. 2(a)). The [OI] 63 µm line is clearly visible while the sensitivity of the instrument allows the detection, at the limit, of the CO (15-14) and CO (17-16) lines (with peaks about 3 - 4\textsigma) but is insufficient to detect the CO (19-18) line. Estimation of the integrated intensities, or their upper limits, can be obtained with reasonable uncertainties (about 24%, see Table 3) since the theoretical FWHM of the lines, related to instrumental characteristics, are known. Emission lines observed by HIFI are well spectrally resolved and with high signal-to-noise ratios (Fig. 2(b)). However, the observing mode chosen for the [CII] 158 µm line has caused a contamination from the background emission of the HII region. The strategy to remove this contamination and extract the [CII] 158 µm line flux is detailed in Appendix A.

In the case of the two Orion proplyds observed with Herschel (HST10 and 244-440), the low-resolution spectroscopy obtained by PACS (Figs. 2(c) and 2(e)) is highly contaminated by the nebula so that the observed emissions mainly come from the background and not from the targets. Only upper limits of the target emission can thus be extracted for the CO (15-14), CO (17-16), CO (19-18) and [OI] 63 µm lines. The high-resolution spectroscopy provided by HIFI enables to spectrally resolve the emission of the Orion proplyds and the emission from the nebula. This is clearly visible with the [CII] 158 µm fine-structure line, on Figs. 2(d) and 2(f), where the targets emission lines are detected at a $v_{\mathrm{LSR}}$ about 21 - 22 km s-1 and located on the wing of the nebula emission line that peaks close to a $v_{\mathrm{LSR}}$ of 10 km s-1. Unfortunately, the velocities of the Orion proplyds shift the [CII] lines at the same frequency of the strongest [13CII] line from the UV-illuminated surface of Orion molecular could (Goicoechea et al. 2015). Here again, only upper limits can thus be extracted. For HST10, the maximum brightness of the [CII] 158 µm line, after correcting the beam dilution, is 4.72 $\times$ 10-6 W m-2 sr-1 which is the same order of magnitude as the one predicted by Störzer & Hollenbach (1999) who probably overestimated it. At low velocities, several emission lines are observed: three components close to a $v_{\mathrm{LSR}}$ of 0 km s-1 and one located about -15 km s-1. Those lines likely correspond to the emission of the 12C+ from the Veil in front of the nebula seen in HI emission and absorption at those velocities (van der Werf et al. 2013). The line at a $v_{\mathrm{LSR}}$ about 70-75 km s-1 is another 13C+ line from the nebula (Ossenkopf et al. 2013). Even with the high spectral resolution of HIFI, no emission from HST10 and 244-440 is visible for the CO (7-6) and CO (10-9) transitions. The peak in the corresponding spectra is the emission from the nebula and the negative signatures are ghosts related to the observing mode (frequency switch, see Table 2). Constraining upper limits of the emission from the targets can still be derived (see Table 3) assuming that the width of the lines should be smaller than the one detected from the 12C+ emission and that the peak of the lines is below the detection threshold, i.e. three times the noise level.

Level-2 calibrated spectra from the PACS line spectroscopy also allow to extract the continuum emission which is underlying the line (see Fig. 2). The continuum extraction was done for the Carina candidate proplyd (105-600) only because the far-infrared signal from the two Orion proplyds is drowned in the emission of the nebula. In a five-by-five spaxels image of PACS, the emission from the target is estimated by extracting the flux in the central spaxel, correcting it for losses to the neighbouring spaxels because of the point spread function, and subtracting the mean observed flux in the 16 external spaxels. The obtained values are presented in Table 4.

2.1.4 Complementary data

For the Carina candidate proplyd 105-600, the dataset is complemented with optical to submillimeter data from other instruments. From the online archives, we retrieved the available data of the Spitzer/IRAC and MIPS (PI: N. Smith, program-ID: 3420, 30848, Smith et al. 2010b) instruments, plus the photometry obtained by Herschel/PACS and SPIRE (PI: S. Molinari, program-ID: OT2_smolinar_7). From this data, we have extracted the continuum emission by direct measurement of the spectral radiance or by aperture photometry (Table 4). We also retrieved the data of the Atacama pathfinder experiment (APEX) 12-m telescope from Sahai et al. (2012) who used the LABOCA and SABOCA bolometers (beam size of respectively 7.8″ and 19″) to measure continuum emission and CHAMP+ and FLASH+ spectrometers (beam sizes from 7.7″ to 18.6″) to observe molecular lines. Finally, we used the optical HST/ACS H\textalpha image (PI: N. Smith, program-ID: 10475, Smith et al. 2010a). Performing photometry at the ionisation front, we obtained a brightness (including H\textalpha and [NII] 6583 Å lines) of $4\times 10^{-7}$ W m-2 sr-1 which has been used to estimate the electron density at the ionisation front of $n_{\mathrm{e}}\leq 680$ cm-3.

For the proplyd 244-440, we also included the OH 84 µm line detected in a PACS survey of the Orion Bar and described in Parikka et al. (2016). HST10 has been observed in the H2 1-0 S(1) line by Chen et al. (1998), and we use their value for the intensity of this line. The [OI] $6300$ Å line was observed with HST and we use the value given in the table 4 of Bally et al. (1998).

2.2 ALMA detection of the 203-506 proplyd

Recently, coauthors of this paper have serendipitously detected intense HCO+ (4-3) emission towards the 203-506 proplyd with ALMA. This detection was quite surprising since only a few, and generally very massive, proplyd silhouette disks have been detected with ALMA in HCO+ and reported in the literature (Williams et al. 2014). Since HCO+ is a PDR tracer, and given that other PDR tracers have been detected towards 203-506 (see below), we have decided to include it in our study. It also complements our sample quite well since it is a silhouette proplyd (Bally et al. 2000), i.e. seen against the bright background in visible. The disk is edge-on and situated just south of the Orion Bar. Its surface is bright in [OI] 6300 Å and a faint [OI] jet is also visible (Bally et al. 2000). It does not exhibit any visible envelope as the others probably because it is lying in the neutral PDR resulting from the photodissociation of the Orion bar. The fact that the ionisation front of the M42 HII region has not yet reached 203-506 suggests that this disk has emerged out of the molecular cloud only very recently and hence that it has been exposed to UV irradiation for a relatively short time compared to classical proplyds. The 203-506 proplyd was detected as part of the ALMA cycle 1 project dedicated to the observation of the Orion bar (PI: J. Goicoechea). The observations were conducted in Band 7 with an angular resolution of 1″ insufficient to spatially resolve the disk. The detailed analysis and data reduction are presented in Goicoechea et al. (2016). 203-506 was clearly detected in the HCO+ (4-3) line. For other tracers with lower critical densities such as low-$J$ CO lines in Band 7 of ALMA, the emission is dominated by the Orion molecular cloud and hence the disk cannot be detected. The derived integrated intensity for the HCO+ (4-3) line is presented in Table 3. Emission from the H2 1-0 S(1) line was also observed with the BEAR instrument (Noel et al. 2005) at the Canada France Hawaii Telescope. This data is presented in Noel (2003) and we use their intensity. The [OI] $6300$ Å line was also observed by HST. The line integrated intensity was extracted from the HST image and corrected from extinction ($A_{\mathrm{V}}=1.26$) using the extinction map of the Orion nebula derived from the MUSE data by Weilbacher et al. (2015). The intensity is reported in Table 3.

3 1D-model of the PDR in an externally-illuminated protoplanetary disk

Since their discovery, proplyds have been extensively studied with the help of PDR models in order to determine their physical properties (e.g. Johnstone et al. 1998; Störzer & Hollenbach 1999; Richling & Yorke 2000; Adams et al. 2004; Clarke 2007; Walsh et al. 2013). This approach is well-suited when trying to explain the FIR lines tracing the warm gas present at the surface layer of the disk and in the envelope. In this section we describe our modelling methodology.

3.1 Proplyds as PDRs

Observations (Sect. 2) and preliminary modelling we performed showed that two different components, with different densities, are necessary to describe a proplyd and correspond very likely to the two distinct regions of the PDR, respectively the disk surface layer and the envelope. Given that, the geometry and parameters that we finally used are summarised in Fig. 3 and Table 5. We consider that a proplyd consists of an envelope (fed by the photoevaporation of the disk) with the apparent shape of an ellipse with semiaxes $a_{\mathrm{env}}$, $b_{\mathrm{env}}$, and defined by a mean density $n_{\mathrm{env}}$, a mean dust temperature $T_{\mathrm{D,env}}$ and $f\mathrm{{}_{C}^{PAH}}$ the fraction of elemental carbon in the gas phase locked in Polycyclic Aromatic Hydrocarbons (PAHs). At the center of the envelope lies the disk with an apparent surface $S_{\mathrm{disk}}$, a mean dust temperature $T_{\mathrm{D,disk}}<T_{\mathrm{D,env}}$ and a mean density $n_{\mathrm{disk}}$ in the surface layer. Note that we are not modelling the properties of the bulk of the disk since we limit ourselves to the PDR which corresponds to the external layers of the proplyd (envelope and disk surface). Finally, the object is irradiated by FUV photons along the major-axis. The corresponding FUV radiation field $G_{\mathrm{0}}$ is expressed in units of the Habing field. These parameters are used to predict both the dust spectral energy distribution (Sect. 3.2), and the gas emission (Sect. 3.3).

3.2 Dust spectral energy distribution

Interstellar dust is mainly composed of small grains (10 nm to 0.1 µm, see e.g. Compiègne et al. 2011) of silicates and carbonaceous compounds. Disks form from this medium and evolve with time in a way that their grains get larger and their dust-to-gas mass ratio increases (see e.g. the review of Williams & Cieza 2011). The dust properties in the photoevaporation flow (envelope) are expected to differ from that of the disk because of a physical selection in size (e.g. Owen et al. 2011). However, we will initially assume that their dust populations are quite primitive and close to the ones of the ISM. We will then discuss the effects of grain evolution (Sect. 3.3 and App. B).

The far-infrared spectrum of proplyds is dominated by the emission of the largest grains from that range, which are at thermal equilibrium with the radiation field. The thermal component of their spectral energy distribution can be modelled knowing the opacity $\tau_{\nu}$ and the mean dust temperature ($T_{\mathrm{D,env}}$, $T_{\mathrm{D,disk}}$444Contrary to the PDR approach which is limited to the surface layer of the disk, the mean dust temperature in the disk $T_{\mathrm{D,disk}}$ is valid for the bulk of the disk.). Indeed, for a given dust temperature $T_{\mathrm{D}}$, the emitted spectral radiance follows the relation

[TABLE]

where $B_{\nu}{}$ is Planck’s law. The dust opacity $\tau{}_{\nu}{}$ at the frequency $\nu{}$, is defined by

[TABLE]

where $\tau{}_{\nu{}_{0}}$ is the dust opacity at the reference frequency $\nu{}_{0}$ and $\beta$ is the spectral index. Following Planck Collaboration et al. (2011), we used a spectral index of 1.8 and a reference wavelength of 250 µm (1200 GHz) for which they found that, in the molecular phase of the ISM, the total column density of hydrogen, $N_{\mathrm{H}}$, is related to the dust opacity by the empirical relation

[TABLE]

The mean column density of hydrogen is related to the total mass of gas and dust by

[TABLE]

where $\Omega{}$ is the observed surface in the instrumental beam knowing the distance to the source, $\mu{}\leavevmode\nobreak\ m_{\mathrm{H}}$ the mean particle mass per hydrogen atom with $m_{\mathrm{H}}$ the mass of one hydrogen atom and $\mu{}=1.4$. Finally, for each component (disk and envelope), the emission can be written as a function of the total mass $M$ and the dust temperature $T_{\mathrm{D}}$,

[TABLE]

The mid-infrared dust emission is dominated by PAHs which are found in the envelope and disk surface of proplyds (Vicente et al. 2013). This component is included in the model using the PAH emission modelled with the DustEM code (Compiègne et al. 2011).We use the standard PAH population given by them, i.e. a log-normal size distribution of neutral molecules and one more for singly-charged molecules. The intensity of PAH emission then only depends on the incident FUV radiation field $G_{\mathrm{0}}$, the abundance of PAHs given as $f\mathrm{{}_{C}^{PAH}}$ and the total mass of the envelope $M_{\mathrm{env}}$.

3.3 Line emission using the Meudon PDR code

To predict emission lines from a PDR, one needs a model that correctly describes geometrical and micro-physical elements. Doing both with accuracy is not currently possible for computational time reasons. In this study, we choose to favor a proper description of the physics and chemistry rather than geometry, also taking into account that we do not have enough angular resolution to perform a careful geometrical analysis. From this perspective, the 1D Meudon code is well suited since it is constantly updated to include latest improvements on the microphysics (e.g. H2 formation and excitation, photoelectric effect, grain temperature fluctuations etc. (Le Petit et al. 2006; Le Bourlot et al. 2012; Bron et al. 2014, 2016). The Meudon PDR code is a 1D-model that considers a cloud made of a stationary plane-parallel slab of dust and gas illuminated by a given radiation field. The model iteratively solves the UV radiative transfer, thermal balance and chemistry at each point of the cloud. Once the computation of atoms and molecules abundances and physical conditions is done for each point in the cloud, the model can calculate the local emissivities of lines. In our case, the integration of this quantity gives us the brightness of the disk surface layer and the brightness of the envelope. Weighted those quantities by the size of the corresponding emitting region enables us to derive the line integrated intensity that can be compared with the spatially unresolved observations that we have.

Here, we use the 1.4.4555http://pdr.obspm.fr/PDRcode.html public version of the code. We use the default values for the dust-to-gas mass ratio, dust properties, cosmic ray ionisation rate and turbulent velocity (Table 6), as well as elemental abundances appropriate for the ISM (see Le Petit et al. 2006). However, grain properties in proplyds are likely to differ from those found in the ISM. Indeed, the smallest grains or PAHs seem to be underabundant in proplyds compared to the ISM: this is observed for HST10 by Vicente et al. (2013) and for 105-600 in this study (see Sect. 4.1.1). This observational fact is in agreement with the low detection rate of PAHs in isolated disks around young low-mass stars (Geers et al. 2006; Oliveira et al. 2010). We set the minimum grain radius to 3 nm in the model, i.e. excluding the smallest particles such as PAHs, for consistency with that observational fact and assuming that this is a general trend for proplyds. In the PDR model, the main impact of this modification is to remove the contribution of these small grains to the photo-electric heating and to the H2 formation rate. This parameter has not a significant impact on the modelled line fluxes because proplyds are dense PDRs where the grain photoelectric effect may not always dominate the heating (this will be discussed in App. B according to our results). Grain growth inside an evolving object tends to modify the grain distribution by increasing the maximum grain radius $a_{\mathrm{max}}$, lowering the absolute value of the power-law index and finally it reduces the FUV extinction curve. Since these parameters can not be constrained accurately based on current observations of proplyds, we have explored their effects by changing independently:

• •

the maximal radius $a_{\mathrm{max}}$ up to 3 µm, consistently with observations of Orion proplyds (Shuping et al. 2003; Vicente 2009a);

• •

the power-law index down to 3.0 (e.g. Testi et al. 2014);

• •

the FUV extinction using various curves of Fitzpatrick & Massa (1988) including a case with a dust FUV extinction cross section per H nucleus $\sigma_{\mathrm{ext}}\simeq 8\times 10^{-22}$ cm2 as used for a sample of proplyds in the Orion nebula by Störzer & Hollenbach (1999).

The effects on the model results are discussed in App. B and the general conclusion is that the changes do not affect the model results significantly. Hence, we keep the default ISM properties throughout this paper.

The input parameters of the Meudon code are the density profile as a function of visual extinction $A_{\mathrm{V}}$, and the incoming FUV radiation field. These parameters are related to those of the proplyd model presented in Sect. 3.1 in the following way. The impinging radiation field is set by the value $G_{0}$ arriving at the surface of the envelope. The density profile between the envelope surface, or ionisation front, at $A_{\mathrm{V}}=0$ and the disk surface at $A_{\mathrm{V}}^{\mathrm{jump}}$ is constant with a value $n=n_{\mathrm{env}}$. At the envelope - disk interface (i.e. at at $A_{\mathrm{V}}^{\mathrm{jump}}$), the density jumps to $n=n_{\mathrm{disk}}$ and is kept constant up to $A_{\mathrm{V}}^{\mathrm{max}}$. Since the parameter $A_{\mathrm{V}}^{\mathrm{jump}}$ gives the position of the envelope - disk interface, it is also the visual extinction of the envelope and is set by

[TABLE]

where $\kappa$ is the parameter which relates the envelope column density $N_{\mathrm{env}}$ to the extinction at the density jump, and can be calculated from Table 6, $\kappa=N_{\mathrm{H}}/A_{\mathrm{V}}=1.87\times 10^{21}$ cm*-2* (assuming ISM properties). The envelope column density is defined by $N_{\mathrm{env}}=n_{\mathrm{env}}\times a_{\mathrm{env}}$. Finally, the value of $A_{\mathrm{V}}^{\mathrm{max}}$ is set to 10. The exact value is not critical but we have carefully checked that it is high enough. Indeed, most of the UV radiation is rapidly absorbed and most of the emission arises from regions of low extinction so increasing this parameter does not impact the results.

3.4 Fitting strategy

In order to estimate the unknown parameters of the model (see Table 5), we are looking for the set of values that best fit the observations based on a chi-square minimisation. For the spectral energy distribution of the dust, the model numerically converges to the best solution in the space of parameters using a Nelder-Meadwith non-linear minimisation script. For the line emission, we run the model hereabove on a grid for $n_{\mathrm{env}}$ and $n_{\mathrm{disk}}$ in a range of realistic values. The apparent size of the disk $S_{\mathrm{disk}}$ is adjusted in the case of 105-600 while is known for the other proplyds. The best model is the one which minimizes the chi-square between observed and predicted line integrated intensities.

4 Results of proplyd models

In this section, we present the results of the comparison between the models and observations. We start with 105-600, for which the largest number of detections is available and follow with the Orion proplyds.

4.1 Proplyd 105-600

4.1.1 Dust spectral energy distribution

The spectral energy distribution of 105-600 is shown in Fig. 4. The results of the fit with the model to the observed SED are given in Table 7. Under the assumptions described in Sect. 3.2, the first component, the disk, is composed of a cold dust population ($\sim$ 20 K) and a total mass (dust and gas) of 0.67 M☉. The envelope corresponds to a hotter dust population ($\sim$ 40 K) and contains a total mass of 0.048 M☉. This mass estimate is compared to other methods in Sect. D.2 which yield similar values both for the disk and envelope. The estimated PAH abundance in the envelope is $f\mathrm{{}^{PAH}_{C}}=0.0024$, about 30 times less than in the PDR of the NGC 7023 reflection nebula (Berné & Tielens 2012). The depletion of PAHs within a proplyd-like object has already been observed in HST10 (Vicente et al. 2013) with a relative abundance three times lower than in 105-600, i.e. about 90 times less than in NGC7023. These low abundances in proplyd-like object could result from PAH destruction by the strong UV radiation field or sedimentation of PAHs in the disk but this is currently still unclear.

4.1.2 Line emission

To compare our models with observations, we have ran a grid of models with density ranges $n_{\mathrm{env}}=1\times 10^{3}-2\times 10^{4}\leavevmode\nobreak\ \mathrm{cm^{-3}}$ and $n_{\mathrm{disk}}=1\times 10^{5}-1\times 10^{7}\leavevmode\nobreak\ \mathrm{cm^{-3}}$, with constant grain properties as defined previously in Table 6 ($a_{\mathrm{min}}=3$ nm, $\delta{}=0.01$, see App. B for the impact of those parameters). Fig. 5 presents the comparison between all the models and the observations. From the \textchi2 analysis (Fig. 5(b)), a minimum is clearly found where the best-fit model is the one with $n_{\mathrm{env}}=8.5\times 10^{3}\leavevmode\nobreak\ \mathrm{cm^{-3}}$, $n_{\mathrm{disk}}=7\times 10^{5}\leavevmode\nobreak\ \mathrm{cm^{-3}}$ and a value $S_{\mathrm{disk}}=2\times 10^{7}$ AU2. Fig. 6 presents the comparison between the observed line integrated intensities and the best-fit model. All lines are reproduced to a factor of a few or better, with the notable exception of HCN (4-3) which is under-predicted by more than an order of magnitude. HCN emission is a tracer of dense gas, whereas our model with a constant density is adopted only for the surface layers of the disk. To go further and improve the modelling of the low-$J$ CO and HCN lines would require a detailed disk model not limited to the PDR as here.

The gas temperature at the edge of the envelope (ionisation front) increases with the envelope density but in a relative small range ($\sim$ 150 to 400 K, see Fig. 5(c)) compared to the variations associated to the grains properties (Fig. 13 and 14). The gas temperature at the surface of the disk mainly depends on the envelope density (see Sect. 5 for more details). The best-fit model corresponds to a value of about 1100 K. Models with slightly lower envelope densities, ans thus slightly higher disk surface temperatures, are still close to the observations while the decrease in disk surface temperature caused by an increase in the envelope density makes models diverge significantly from observations. This is strong evidence that the disk surface is hot with temperature around 1000 K or more.

4.1.3 Physical structure of the PDR

As an output of the model it is possible to obtain the 1D physical structure of the PDR in the studied proplyds. Note that this structure is only limited to the PDR, i.e. the envelope and the surface layer of the disk, and so that the extracted physical properties are not valid deeper in the disk where much larger densities are expected and UV photons are strongly attenuated. The 1D profiles for the best model of 105-600 are shown in Fig. 7. It can be seen that the H/H2 transition, defined here as the region where 10 to 90% of H is in H2, encompasses the envelope - disk interface ($A_{\mathrm{V}}^{\mathrm{jump}}=0.742$, see eq. (6)) ranging from $A_{\mathrm{V}}=0.729$ (inside the envelope, right before the interface) to $A_{\mathrm{V}}=0.748$ (inside the disk). The C+/CO transition starts inside the molecular disk, at $A_{\mathrm{V}}=0.742$ and stops at $A_{\mathrm{V}}=2.49$.

The gas temperature is 285 K at the ionisation front and increases up to 650 K at $A_{\mathrm{V}}=0.63$ (Fig. 7, second pannel). It then decreases to about 320 K before jumping to 1100 K at the disk surface and then slowly decreases. Note that a smoother density profile, rather than a jump, does not significantly change the temperature profile except that it is slightly smoothed too. We thus choose to keep a jump for simplicity. The initial increase inside the envelope can be attributed to an increase of the photoelectric-effect efficiency related to the increasing concentration of neutral grains (Fig. 7, third pannel). The trend is inverted just before reaching the disk (and the H/H2 transition) where the efficiency drops. The photoelectric-effect is by far the main heating mechanism inside the envelope. However, at the disk surface, several heating mechanisms contribute significantly (Fig. 7, third pannel): chemical reactions (37%), photoelectric-effect (24%), formation of H2 on grains (17%,) collisional de-excitation of UV-pumped H2 (15%) and direct ionisation/dissociation by FUV photons (7%). Beyond $A_{\mathrm{V}}=1$, the photoelectric-effect becomes again the main heating mechanism up to $A_{\mathrm{V}}=4.5$, where the gas-grains collisions take the lead and keep the gas temperature close to the dust temperature. However, at this location, our simple constant density PDR model is not appropriate anymore and proper modelling will require specific disk models.

Inside the envelope, owing to the relatively high gas density, the main cooling mechanism is the fine-structure line emission of [OI] 63 µm, which is one order of magnitude stronger than the [CII] 158 µm line. At the surface of the disk, the cooling is shared between the cascade of H2, the gas-grains collisions, the emission of O and CO. Deeper in the disk, CO and O dominate the cooling.

4.2 Results for Orion proplyds

For the three Orion proplyds, few detections of gas lines are available. For these sources the far-infrared dust emission was not detected since it is drowned in the emission from the nebula. We therefore only attempt to reproduce the observed line emission with our model. Intense emission of the [OI] 6300 Å line is observed towards the disks of HST10 and 203-506 as reported by Bally et al. (1998). This line arises from the hot neutral layers of the PDR. Its excitation is the result of the rapid formation of OH through the reaction of O with H2, followed by photodissociation leaving a large fraction of O in an electronically excited state which radiatively decays through emission in the 6300 Å line. Since this process was not included in the public version of the Meudon PDR code (it will be in the next one), we therefore calculate the intensity of the [OI] 6300 Å line using the prescription of Störzer & Hollenbach (1998), using the column density of OH predicted by our model. In order to fit the observations, we follow the same strategy as for 105-600, i.e. we vary the envelope and disk density. Comparison between observations and the best-fit model are shown in Fig. 8, while parameters and outputs of the best-fit models are given in Tables 5 and 8.

For HST10, models cannot produce more than $10^{-6}$ W m-2 sr-1 of [CII] emission – including changing the envelope density, the small grains abundance or dust extinction properties – which is well below the upper limit extracted from observation. This is consistent with the idea that the observed emission should be dominated by the [13CII] emission from the HII region expected at the same frequency. Overall this suggests that the envelope properties are not well constrained here. The [OI] 6300 Å line intensity is initially strongly underestimated but is better reproduced if the mechanism of formation for the excited oxygen described by Störzer & Hollenbach (1998) is included in the model and if an $n_{\mathrm{env}}\leq 4\times 10^{5}$ cm-3 is used. Above this value, the envelope column density is large enough that OH starts to form inside the envelope, which is incompatible with observations showing that the emission occurs close to the disk surface (Bally et al. 1998). An envelope density of $n_{\mathrm{env}}\leq 4\times 10^{5}$ cm-3 corresponds to an envelope column density of about $N_{\mathrm{env}}=3\times 10^{21}$ cm-2 which is slightly below the estimated value of $5.5\times 10^{21}$ cm-2 from Störzer & Hollenbach (1998). The H2 and [OI] 6300 Å lines are best explained with a disk density of the order of $10^{6}$ cm-3. The best-fit model gives $n_{\mathrm{env}}=2.5\times 10^{5}$ cm-3 and $n_{\mathrm{disk}}=4\times 10^{6}$ cm-3. The disk density is consistent with the estimate of Störzer & Hollenbach (1998) who found the same value. The temperature at the disk surface predicted by the best fit model is 1150 K.

For 244-440, [CII] arising mainly from the envelope is probably detected but, as HST10, drowned in the nebula emission. The OH line is well detected and should arise from the the disk (under the assumptions that this line is being emitted by the proplyd and no other object of the Orion Bar in the line of sight). With two lines, whose one is an upper limit, the two densities will be poorly constrained, but we nevertheless found a clear best-fit for a model with $n_{\mathrm{env}}=1\times 10^{5}$ cm-3 and $n_{\mathrm{disk}}=2\times 10^{6}$ cm-3. The temperature at the disk surface predicted by the best fit model is 1270 K. These values are consistent with those found for the other Orion proplyds.

For 203-506, the disk density is well constrained by the observation of the H2 and HCO+ lines and the best-fit value is $n_{\mathrm{disk}}=1\times 10^{6}$ cm-3. The ”envelope” density that corresponds to the ambiant PDR resulting from the photoevaporation of the Orion bar (Goicoechea et al. 2016) is estimated to be $n_{\mathrm{env}}=4.5\times 10^{4}$ cm-3. The temperature at the disk surface (which we consider being at the H/H2 transition) predicted by the best fit model is 1080 K. Here again, we need to include the mechanism of Störzer & Hollenbach (1998) in the model in order to get closer to the observed [OI] 6300 Å line intensity but is still highly underestimated. This discrepancy could be explained by the emission from shocked regions or that another mechanism to form OH is missing, e.g. by reactions with vibrationally excited H2.

The physical structures extracted for the Orion proplyds are not shown here because they are less well constrained than 105-600. They are found similar to the one given for 105-600 in Fig. 7, with a rapid increase of the gas temperature at the disk surface and the H/H2 transition layer located close to it.

4.3 Comparison to other models

It is interesting to compare the results obtained here with those of earlier studies of proplyds using PDR models. When doing this comparison it is important to keep in mind that these studies usually do not compare their predictions to actual observations of PDR tracers towards these objects.

Störzer & Hollenbach (1999) and Adams et al. (2004) use a PDR model based on the classical Tielens & Hollenbach (1985) model. Contrary to what is presented here, they consider a single constant density structure. Störzer & Hollenbach (1999) compare their models to the emission of vibrationally excited H2 and [OI] 6300 Å emission observed in Orion proplyds (including HST10) and find good agreement.They also predicted the far infrared line intensities of the [CII] line, for which they give values ranging between 2.4 and $3.1\times 10^{-6}$ W m-2 sr-1, which assuming an emitting size equal to the size of the envelope projected on the plane of the sky ($2.6\arcsec\times 1\arcsec$, see Table 1) corresponds to a flux of 1.15 - 1.50 $\times 10^{-16}$ W.m-2, just below the high upper limit of $2.3\times 10^{-16}$ W.m-2 that we obtained. While not excluded, their value is probably overestimated since the observation is supposed to be highly contaminated by the nebula emission. They have also predicted strong [OI] 63 µm emission close to our high upper limit of detection leading to the same possible conclusion. Detailed observation are needed to conclude clearly. The study of Adams et al. (2004) is essentially theoretical. While focusing on different FUV field ranges, the temperature profile is characterized in both studies by a region of constant and high temperature (¿ 4000 K) at low column densities followed by a drop of temperature at the location of the H/H2 transition, which is situated, in their models, just above the disk. The temperature profile which we find here (Fig. 7) differs greatly from that of Störzer & Hollenbach (1999) and Adams et al. (2004) in that the temperature we find at low column densities is much lower (typically a few 100 K). This is mainly due to two fundamental differences which reduce the heating in the envelope in our models: 1) the lower density does not allow collisional de-excitation of H2 to be important, 2) the low abundance of small grains that we have considered (in agreement with the observation that PAHs are underabundant in HST10 and in 105-600 (Vicente et al. 2013, and Sect. 4.1.1), reduces photoelectric heating drastically. Nevertheless, an important aspect in the frame of photoevaporation is that the temperatures at the disk surface derived by Störzer & Hollenbach (1999) and Adams et al. (2004) are typically of the order of a few 100 K, i.e. lower than ours. Our prediction for the H2 (1-0) S(1) line for HST10 is in good agreement with theirs and with our observations.

5 Effect of the H/H2 transition layer location on the gas thermal balance at the disk surface

In this section, we explore the effect of the location of the H/H2 transition layer on the gas thermal balance, and the gas temperature, at the disk surface, i.e. the density jump in our model. As it is described in the book of Tielens (2005), the location of this layer, i.e. the column density of hydrogen nucleus, is proportional to the FUV flux $\mathrm{G}_{0}$ and inversely proportional to the density. Since the column density between the ionisation front and this layer is mainly set by the envelope density in the case of proplyds, we run our model on a grid of values for $n_{\mathrm{env}}$ and $\mathrm{G}_{0}$ from $2\times 10^{3}$ to $2\times 10^{5}$. We also considered different envelope sizes: a large one corresponding to the size of the biggest studied proplyd, 105-600, and a small one corresponding to the size of HST10.

5.1 Temperature

First of all, we extracted the gas temperature at the disk surface for each model of the grid to see its evolution with the position of the H/H2 transition layer. This is illustrated in Fig. 9 which presents the gas temperature at the disk surface as a function of the position of the H/H2 transition which is captured by the ratio $P_{\mathrm{trans}}=A_{\mathrm{V}}^{\mathrm{jump}}/A_{\mathrm{V}}^{\mathrm{H-H_{2}}}=N_{\mathrm{env}}/N_{\mathrm{H-H_{2}}}$, where $A_{\mathrm{V}}^{\mathrm{H-H_{2}}}$ is the visual extinction, or position, at which the transition layer starts, i.e. where 10% of the hydrogen is in molecular form, and $A_{\mathrm{V}}^{\mathrm{jump}}$ is the visual extinction setting the position of the density jump corresponding to the disk surface in our model. The H/H2 transition is thus totally located inside the disk when $P_{\mathrm{trans}}<1$, and is located or starts in the envelope when $P_{\mathrm{trans}}>1$. A similar behaviour is observed for all the studied cases: a hot ($T=1500$ - $5000$ K) regime prevails when the transition is located inside the disk ($P_{\mathrm{trans}}<1$) while a warm regime ($T<500$ K) is observed when the transition is located inside the envelope ($P_{\mathrm{trans}}>1$). A transitional regime ($T\sim 500$ - $1500$ K) lies between the hot and warm regimes at $P_{\mathrm{trans}}\sim 1$.

5.2 Heating and cooling mechanisms

To gain insight into the mechanisms driving this temperature behaviour, we also studied the evolution of the heating and cooling mechanisms as a function of the H/H2 transition layer position. Results are given in Fig. 10 and Fig. 11 respectively. The following observations can be made: when the envelope density (and visual extinction) is small and the transition is inside the disk ($P_{\mathrm{trans}}<1$, $T\sim 1500$ - $5000$ K), the collisional de-excitation of UV-pumped H2 is very efficient and dominates the heating, except for very intense FUV fields ($G_{\mathrm{0}}>10^{5}$) where the heating by the photo-electric effect can be even higher, while cooling is dominated by gas-grain collisions. This corresponds to the hot regime. When the transition is located inside the envelope ($P_{\mathrm{trans}}>1$, $T<500$ K), the photoelectric effect on grains dominates the heating, and cooling is dominated by emission in fine structure lines of oxygen. This corresponds to the warm regime which is the classical solution for thermal balance usually observed in PDRs. Finally, in the transitional regime, i.e. when the transition is close to the surface but starts inside the envelope ($P_{\mathrm{trans}}\sim 1$, $T=500$ - $1500$ K), the contribution of a collection of heating processes (exothermal reactions, H2 formation and direct absorption of ionizing/dissociation photons, exothermal chemical reactions) control the heating and cooling also occurs through multiple processes (e.g. collisional excitation of H2, OH and CO emission lines)

5.3 Evidence for a specific energetic regime

At this point it is interesting to note that for each object, the H/H2 transition in the best-fit model is located close to the disk surface ($P_{\mathrm{trans}}\gtrsim 1$, cf. Table 8 and Fig. 9). This is in agreement with spatially resolved observations which for instance show that the emission of vibrationally excited H2 (for HST10) and [OI] 6300 Å (for HST10 and 203-506) arises from the disk surface or just above (Chen et al. 1998; Bally et al. 1998). This implies that proplyds appear to be in the transitional regime with temperatures at the disk surface close to 1000 K. Guided by Fig. 9, we build the following reasoning: when the isolated disk (before the formation of the envelope) is irradiated by the external FUV photons, the H/H2 transition forms inside the disk, and the temperature rises to some several 1000 K, driving an intense photoevaporation flow. This increases the column density of the envelope $N_{\mathrm{env}}$, dragging the H/H2 transition out of the disk ($P_{\mathrm{trans}}=N_{\mathrm{env}}/N_{\mathrm{H-H_{2}}}>1$) which lowers the temperature at the disk surface. Consequently, the evaporating mass flux will decrease and so will the envelope column density due to spherical divergence. At one point the H/H2 transition will go back inside the disk ($P_{\mathrm{trans}}<1$), temperature will rise, mass-loss increase as well as envelope column density, etc.

Overall, this suggests that proplyds are in an equilibrium where the temperature of the disk surface is set so that mass loss is sufficient to keep the envelope column density to a value that maintains the H/H2 transition close to the disk surface. In this specific case, the disk surface temperature is set by a combination of numerous heating and cooling processes and reaches values of the order of 1000 K (for radiation fields between $G_{\mathrm{0}}$=2000 and $G_{\mathrm{0}}$=$2\times 10^{5}$).

6 Photoevaporation and mass-loss rates

At this point it is interesting to discuss photoevaporation and mass loss of proplyds in the context of our results. Namely, knowing the size of the disks and a possible range of masses for the central stars from independent studies, we can use the surface temperature and densities from our models (Tables 5 and 8) to estimate photoevaporation rates.

6.1 Definitions

The critical radius where the kinetic energy of gas particles is equal to the gravitational potential is defined as (Hollenbach et al. 1994)

[TABLE]

where $c_{\mathrm{S}}$ is the sound speed, i.e.

[TABLE]

The term $\mu\leavevmode\nobreak\ m_{\mathrm{H}}$ is the mean particle weight, $\gamma$ is the adiabatic index taken as 7/5 here - the value for a perfect diatomic gas - and $\mathrm{k_{B}}$ the Boltzmann constant. The mass loss rate is defined for two cases, which depend on the value of the disk radius $r_{\mathrm{disk}}$:

Supercritical case when $r_{\mathrm{disk}}>r_{\mathrm{g}}$.

The mass loss rate occurs mostly from the disk surface (vertical mass-loss) and is given by (Johnstone et al. 1998),

[TABLE]

Subcritical case when $r_{\mathrm{disk}}<r_{\mathrm{g}}$.

The mass loss rate occurs from the disk edge (radial mass-loss) and is given by (Adams et al. 2004),

[TABLE]

where $C_{0}$ is a dimensionless constant of order unity and $N_{C}$ is defined as the atomic region of the PDR in their models which must be of the order of $N_{\mathrm{env}}$ in ours.

6.2 Carina candidate proplyd 105-600

For 105-600, the disk (or candidate disk) size is unknown. The value of $S_{\mathrm{disk}}$ derived from the model can be converted in an equivalent radius considering that the disk is in the plane of the sky. This yields a value of about 1.1″or 2500 AU. The temperature at the surface of the disk is estimated to be about 1100 K. For such a large radius, $r_{\mathrm{disk}}>r_{\mathrm{g}}$ for masses of central star up to 30 M⊙, and the supercritical regime prevails everywhere. We derive for this case a mass-loss rate $\dot{M}=7.5\times 10^{\mathrm{-5}}$ M${}_{\sun}/$year. The mass-loss at the ionisation front can be derived from the H\textalpha emission, and we find a value of the order of $2\times 10^{\mathrm{-6}}$ M${}_{\sun}/$year, about 40 times smaller. This discrepancy could be explained if the envelope is in a phase where it is gaining mass (see Sect. 5.1 ) while depleting the disk. We note that the Orion proplyds, which seem to be at equilibrium (mass-loss at the ionisation is roughly equal to the ”disk” mass loss computed here, see Sect. 6.3), generally exhibit a tear-drop shaped envelope as HST10 (Fig. 1(b)) as the consequence of this equilibrium. It is thus possible that the sharper shape of 105-600 (Fig. 1(a)) results from an envelope that is still growing in mass. Since the growth phase is expected to be quick, e.g. forming a 0.05 M☉ envelope with such a rate would take few hundred years, 105-600 would thus be a very young object. Alternatively, it is possible that the disk radius we consider here is too large by about an order of magnitude. This would however pose a serious problem to our model when trying to reproduce the high-$J$ CO lines. Since a protoplanetary disk with a radius of 2500 AU appears unrealistic, what we call disk in this case could most likely refer to a young forming disk from the original cloud material (similar to other large ”pseudo-disks” observed, e.g. as in Quanz et al. 2010). Further observations at (sub)millimeter wavelengths to resolve the structure and kinematics of the disk will be required to settle this issue.

6.3 Orion Proplyds

For HST10 and 244-440, Henney & O’Dell (1999) derived the mass-loss rates at the ionisation front from spectroscopic observations of ionized gas tracers and found values of $\left(7\pm 5\right)\times 10^{-7}$ and $\left(15\pm 7\right)\times 10^{-7}$ M☉ yr-1 respectively. For 203-506, since the envelope ionisation front is not observed because the disk is embedded in neutral gas, the mass-loss cannot be derived from ionised gas tracers.

For the Orion proplyds whose disk outer radii are known (Vicente & Alves 2005), we can derive the mass-loss rates from our model using the analytical prescription presented above. The masses of their central stars are not precisely known but taking into account the stellar mass function of the cluster (e.g. Hillenbrand & Carpenter 2000), they are very likely about 0.1 M☉ or less. A surface temperature of the order of 1000 K, as obtained from our study, would result in a gravitational radius of 20 AU, or less, and implies the supercritical case ($r_{\mathrm{disk}}>r_{\mathrm{g}}$) for all disks. The subcritical case is still possible, but unlikely because we need to assume that central stars have a mass of about 1 M☉, resulting in a gravitational radius of the order of 200 AU. We nevertheless consider both cases in our calculations. The mass-loss rates for each case are given in Table 8. In the subcritical regime, we find values of the order a few $10^{-8}$ M*☉/yr, in agreement with the results of Adams et al. (2004), but which are too low compared to the observed values. On the other hand, the supercritical values of mass-loss are in better agreement with the mass loss derived at the ionisation front. Hence, to reconcile the photoevaporation rates derived from our best fit models with those derived from the observations by Henney & O’Dell (1999), one must assume that the flows are in a supercritical regime. This is compatible with the idea that central stars have masses well below one solar mass. If this is indeed the case, given that disk masses do not significantly exceed $0.01$ M☉* in mass, this implies short lifetimes for these disks (of the order of several 104 years) which suggests that we should not see so many of them, or that the ionizing star \textTheta1 Ori C is particularly young, a conundum that is known as the ”proplyd lifetime problem” (Henney & O’Dell 1999). Clarke (2007) studied in details the dynamical evolution of proplyds assuming they are in the subcritical regime and found that their lifetime could in this case be relatively large, solving this issue. In the light of the results presented here it appears that the subcritical assumption may not be correct. Instead our results are in line with the idea that the Orion Nebula is young, i.e. that $\Theta^{1}$ Ori C is less than a few 104 years old, and that proplyds are indeed rare objects which we happen to see in Orion by chance. Another possibility is that the supercritical regime is not sustained at the same rate during the whole proplyd life and declines with time so that the lifetimes could be not so short.

7 Conclusions

The main results presented in this paper are:

• •

We have presented the first far-infrared observations of dense PDR tracers emitted by proplyds, obtained with the Herschel space observatory and ALMA.

• •

Based on the detailed Meudon PDR code, we have developed a 1D-model for the PDR of a proplyd to predict the line emission arising from it.

• •

This model successfully reproduces most lines, with the exception of HCN (4-3) in 105-600. The model is also somewhat lower than the observed intensities of [OI] 6300 Å (and possibly [CII] 158 µm) in HST10 and 244-440. This tracer is most sensitive to the densest gas (or the most diffuse for [CII] 158 µm) in proplyds which are not well described in our simplified model.

• •

For all the sources, at the disk surface, we find densities of $0.7$ to $4\times 10^{6}$ cm*-3* and temperatures of 1000-1300 K.

• •

We find that the position of the H/H2 transition layer is the critical parameter determining the disk surface temperature.

• •

Our results suggest that proplyds are in a self-regulated regime where the temperature of the disk surface is set so that mass loss keeps the envelope column density sufficient to bring the H/H2 transition out of the disk, yet close to its surface.

• •

We find that gas energetics in this specific case results from a complex combination of several heating and cooling mechanisms, hence it can not be captured by the classical solution where photoelectric effect on grains dominates the heating and [OI] and [CII] far infrared emission dominates the cooling.

• •

We derive mass loss rates that are large (of the order of a few 10*-7* M*☉/yr for the Orion proplyds and a few 10-5* M*☉*/yr for 105-600). Most of them are consistent with spectroscopic studies of the ionisation front in theses objects and suggest that proplyds undergo supercritical photoevaporation and indeed may have a short lifetime.

• •

To understand proplyd evolution, in a next step, it is necessary to couple hydrodynamic simulations of a supercritical evaporation flow with a correct treatment of PDR physics and chemistry.

• •

Concomitant with the development of models, spatially resolved observations of molecular lines close to the surface of the disks would help to constrain the photoevaporative flows and to understand the evolution of a proplyd. Especially, the high angular resolution provided by ALMA and the future JWST will enable us to resolve photoevaporative flows and study their properties and dynamics.

Acknowledgements.

This work was supported by the CNRS program ”Physique et Chimie du Milieu Interstellaire” (PCMI). The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) ERC-2013-SyG, Grant Agreement n. 610256 NANOCOSMOS. Javier R. Goicoechea also thanks the Spanish MINECO for funding support under grant AYA2012-32032. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00352.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Appendix A 105-600: case of the [CII] emission

As can be seen on Fig. 2(b), the fine-structure line of [CII] from the target does not seem well defined. Indeed, since the observing mode chosen was dual beam switch (see Table 2), the final spectrum is the result of a spectra obtained on the target (ON position) and others off the target (OFF position). Rather being flat, The off position is apparently not free of emission for this line so that we sum a positive line with a negative one. Fortunately, we have some information about the target as well as for the nebula which is seen in the OFF position. It is known from Sahai et al. (2012) and from our other resolved lines (see Table 3) for this object that it has a $v_{\mathrm{LSR}}$ about $-\leavevmode\nobreak\ 23$ km s-1. This confirm that the positive peak seen on the resulting spectrum is the target emission. Concerning the OFF emission from the nebula, it is known that it comes from a HII region where the temperature is about 104 K. At this high temperature, the width of the line is dominated by the thermal doppler broadening. With the constraints on the position of the target line and the width of the nebular line, we can fit the spectrum (Fig. 12).

Our model with these two gaussian lines fit quite well the part of the spectrum with $v_{\mathrm{LSR}}<-10$ km s-1 while there are still some features on the remaining part which do not overlap the target emission. The constraints are sufficient to retrieved the target profile with a good uncertainty (see Table 3).

Appendix B Effect of grain parameters (test on 105-600)

While grains represent only a small fraction of the total mass, a correct treatment of their role is crucial since the photoelectric effect on dust grains (Bakes & Tielens 1994) is supposed to be one of the major sources of heating and because they are responsible for the UV extinction in the PDR. The grain population is generally not precisely constrained in an evolving protoplanetary disk. It should slightly differs from the one of the ISM, so we investigated the effect of changing the main grain parameters in the model: the grain size distribution given by the minimum radius ($a_{\mathrm{min}}$), the maximum radius ($a_{\mathrm{max}}$) and the power-law index, and the dust-to-gas mass ratio $\delta{}$.

If, contrary to the models presented in the main text, we include the smallest particles, i.e. we set $a_{\mathrm{min}}$ to 0.4 nm instead of 3 nm, we observe very minor changes in the modeled line intensities: low-$J$ CO lines are slightly lower and the [OI] 63 µm is slightly higher (Fig. 13). In the absence of strong dependency, those lines cannot significantly constrain with confidence this grain parameter. The temperature at the surface of the disk is almost constant but the thermal profile is significantly changed in the envelope with strong increase of the temperature when small grains are included because they improve heating by photoelectric effect. According to the electron density estimated at the ionisation front, $n_{\mathrm{e}}\leq 680\leavevmode\nobreak\ \mathrm{cm^{-3}}$, and assuming a pressure equilibrium (constant $P=nT$) between the neutral and ionised gas, we can derive the temperature on the side of the envelope. Note that, rigorously, the pressure in the envelope is likely somewhat larger than in the ionised gas in order to maintain the pressure gradient of the photoevaporation flow. However, the temperature on the side of the envelope should be

[TABLE]

where $T_{\mathrm{HII}}=10^{4}$ K is the temperature in the ionized region. In the case of the best-fit model for 105-600 where $n_{\mathrm{env}}=8.5\times 10^{3}\leavevmode\nobreak\ \mathrm{cm^{-3}}$, the temperature is thus estimated to be $T_{\mathrm{env,IF}}\leq 800$ K. Models are consistent with this result only if small grains are removed by setting the minimal radius to 3 nm or a bit less, i.e. without PAHs. This is consistent with the observed low abundance of PAHs (Sect. 4.1.1).

We investigate the effect of changing the dust-to-gas mass ratio by varying it up to an order of magnitude (Fig. 14). In that case, line fluxes are very sensitive and the best model is the one with the standard value of 100 which is, moreover, compatible with the mass estimates (see Sect. D.2). Assuming a slightly higher dust-to-gas mass ratio is still possible but it will increase the temperature at the ionisation front and will not be consistent with the 800 K limit for the envelope side. We therefore adopt gas-to-dust mass ratio of 100 ($\delta=0.01$) and $a_{\mathrm{min}}=3$ nm.

We also investigated the effect of the grain growth by increasing the maximum radius and decreasing the power-law index.Adding large grains does not significantly impact the photo-electric effect but the FUV extinction curve is changed. To better describe this effect we ran tests using a standard ISM size distribution ($a_{\mathrm{max}}=0.3$ µm and a power-law index of 3.50), but with lower extinctions with the curves measured towards HD38087, HD36982, and HD37023. These span a wide range of dust FUV extinction cross section per H nucleus $\sigma_{\mathrm{ext}}$ including a value as low as $\simeq 8\times 10^{-22}\leavevmode\nobreak\ \mathrm{cm^{-2}}$ which is the value observed towards \textTheta1 Ori C and appropriate for Orion proplyds (Störzer & Hollenbach 1999). The main effect of lowering the FUV extinction is to shift the H/H2 transition deeper in the PDR. This does not change the main heating/cooling processes involved nor the general results of our study. However, this could result in a slightly higher mean density and subsequent total column density for the envelope (as a comparison, for 105-600, the best-fitted density is $n_{\mathrm{env}}=8.5\times 10^{3}$ cm-3 for the Galaxy curve and $n_{\mathrm{env}}=9.0\times 10^{3}$ cm-3 for the one of HD37023). Since, we did not obtain better fits to our observations by changing the extinction curve, we kept a standard Milky Way extinction curve.

Appendix C 105-600: the possibility of emissions from shocks

The Carina candidate proplyd exhibits a bipolar collimated jet that is visible in HST/ACS H\textalpha image (Smith et al. 2010a). Shocks along jets around protostars at comparable evolutionary stages are known to be sources of optical and FIR emission lines (see e.g. Podio et al. 2012). We should thus verify if some of the observed lines could originate in shocks rather than in the PDR.

Fast and dissociative J-type shock models (velocity between 30 and 150 km s-1) of Hollenbach & McKee (1989) could explain a fraction of the observed atomic and molecular lines but predict, for a pre-shock density of $10^{4}$ cm-3 or more, that the shocked surface should be bright in the H\textalpha image which is not the case (see Fig. 1(a)). For lower densities, the emissions in the FIR become negligible compared to the observations. We conclude that J-shocks are unlikely to be responsible for the observed infrared lines. Moreover, except if the jet is almost perfectly in the plane of the sky, the relatively low FWHM of the [CII] line (6.4 km.-1) suggests that this line is not coming from a high-velocity shock.

For low-velocity shocks ($\leq 30$ km s-1) associated with a jet, we looked at the estimated line fluxes from C and J-shock models of Flower & Pineau Des Forêts (2010). C-shocks produce negligeable [OI] emission compared to the observed one while J-shocks are not able to produce the observed [OI] emission without being, at least, two orders of magnitude brighter in CO compared to the observed emission. Atomic lines are thus very unlikely produced by such shocks. Finally, we conclude that atomic lines arise from the PDR of 105-600.

The very low FWHM of CO lines observed with Herschel/HIFI or APEX rule out the origin of this emission from relatively high-velocity shocks associated with the jet. The presence of the jet, however, suggest an accretion of matter onto a young star (see e.g. Konigl & Pudritz 2000). Different low-velocity shocks can be expected in forming disk such as an accretion shock at the surface of the molecular disk, shock waves launched deeper in the disk or residual shocks from the interaction of the jet surface with the disk (Boss & Durisen 2005). We used a C-shock code assuming a shock velocity of 4 - 6 km s*-1*, a magnetic field of 1 - 2 mG and a pre-shock density of $10^{6}\leavevmode\nobreak\ \mathrm{cm^{-3}}$ to estimate the radiated energy to compare it with the UV energy absorbed. In the most favorable case, the radiated energy flux, or radiated power per unit of shocked surface, reaches $8.65\times 10^{-5}\leavevmode\nobreak\ \mathrm{J\leavevmode\nobreak\ m^{-2}\leavevmode\nobreak\ s^{-1}}$. The photo-electric effect injects a fraction $\epsilon{}_{\mathrm{pe}}\approx 1\%$ of the UV flux ($G_{\mathrm{0,mol}}$ times the Habing’s field) into the gas. Hence, the power per unit of disk surface is $\epsilon{}_{\mathrm{pe}}\leavevmode\nobreak\ G_{\mathrm{0,mol}}\leavevmode\nobreak\ 1.6\times 10^{-6}\leavevmode\nobreak\ \mathrm{J\leavevmode\nobreak\ m^{-2}\leavevmode\nobreak\ s^{-1}}$ where $G_{\mathrm{0,mol}}$ is the UV radiation field incoming at the surface of the disk which is about 10% of the initial flux incoming on the envelope for the range of studied envelope densities. With these values and $G_{\mathrm{0,mol}}=2.2\times 10^{3}$, the energy flux is $3.52\times 10^{-5}\leavevmode\nobreak\ \mathrm{J\leavevmode\nobreak\ m^{-2}\leavevmode\nobreak\ s^{-1}}$. To have an important contribution from the shock, the comparison of the two sources of energy implies that the shocked surface has to be about 45% of the disk surface. Moreover, if it happens, some of the emitted lines can be of the order of the observations only if the shock is not FUV-illuminated. Otherwise the emission drops significantly. Hence we cannot fully exclude that shocks could have a contribution to the CO emission, however this requires that the shocked surface is large and that it is protected from UV radiation. Both constraints seem difficult to concile. Altogether, we conclude that shocks are unlikely to play a major role in the observed emission.

Appendix D 105-600: mass and mass-loss rate

D.1 Molecular mass from 12CO and 13CO (3-2) lines

Column density of hydrogen in the molecular disk, and so its mass, can be estimated from the 12CO (3-2) and 13CO (3-2) transition lines. Under the assumption that the 13CO line to be optically thin, the column density of the upper energy level (here 3), noted $N_{\mathrm{u}}$, is linked to the integrated intensity of the line,

[TABLE]

where $\nu$ is the frequency of the transition, $A_{\mathrm{ul}}$ the Einstein emission-coefficient between the upper and the lower levels of the transition. $h$ is the Planck’s constant, $c$ the vacuum speed of light, $k_{\mathrm{B}}$ the Boltzmann’s constant and $\int{T_{\mathrm{B}}\left({}^{13}\mathrm{CO}\right)d\nu}$ the integrated intensity of the line. Additionally assuming that Local Thermodynamical Equilibrium (LTE) conditions are verified for 13CO, the column density of the upper level, $N_{\mathrm{u}}$, is linked to the total one, $N$, by

[TABLE]

where $g_{\mathrm{u}}$ and $E_{\mathrm{u}}$ are respectively the statistical weight and energy of the upper level. $Z$ is the partition function and $T_{\mathrm{ex}}$ is the excitation temperature of the transition. This latter can be calculated thanks to the optically thick emission of 12CO since we have

[TABLE]

where $T_{\mathrm{bb}}=2.7$ K is the background temperature.

From Sahai et al. (2012), we have $T_{\mathrm{B}}\left({}^{12}\mathrm{CO}\right)=33.0$ K so we get $T_{\mathrm{ex}}=40.8$ K. With a line integrated intensity of $\int{T_{\mathrm{B}}\left({}^{13}\mathrm{CO}\right)d\nu}=0.96$ K.km s-1, we find $N\left({}^{13}\mathrm{CO}\right)=3.53\times 10^{14}$ cm-2. Supposing an abundance ratio $\left[\mathrm{H_{2}}\right]/\left[{}^{13}\mathrm{CO}\right]=7\times 10^{5}$, we obtain the column density of $H\mathrm{2}$, $N_{\mathrm{H_{2}}}=9.84\times 10^{22}$ cm-2. The molecular mass, $M_{\mathrm{mol}}$, is finally calculated by

[TABLE]

where $\mu=2.8$ is the mean molecular weigh, $m_{\mathrm{H}}$ is the mass of a hydrogen atom and $S$ is the apparent surface visible is the beam of the instrument at the distance of the observed object. In our case, the molecular mass is $M_{\mathrm{mol}}\approx 0.187\leavevmode\nobreak\ \mathrm{M_{\sun}}$.

D.2 Consistency of mass estimates

We can roughly estimate of the disk mass assuming spherical geometry and making the strong hypothesis that the density is constant throughout the disk. In reality, the density deep in the disk could be much higher so that this estimate can be seen as a lower limit. From our PDR best-fit model, this mass is $M_{\mathrm{disk}}\gtrsim 0.19$ M☉. In the same way, we can evaluate the mass of the envelope assuming that the 3D shape is the one of an ellipsoid with axis of $9.5\arcsec\times 3.7\arcsec\times 3.7\arcsec$ (from Table 1) which gives $M_{\mathrm{env}}\approx 0.028$ M☉ for the best-fit model.

We have used different methods to estimate the mass of each component in the Carina candidate proplyd. The mass of the atomic envelope is estimated to be about $0.05$ M☉ from the SED (Sect. 4.1.1) and about $0.03$ M☉ from the PDR models (Sect. 4.1.2). For the mass of the disk, the SED gives a value of about $0.7$ M☉ (Sect. 4.1.1), consistent with the lower limit of $0.19$ M☉ from the PDR model (Sect. 4.1.2). The fit on the thermal components of the dust emission using modified Planck’s functions is based on empirical parameters such as the spectral index $\beta$ and the efficiency section $\sigma$. The latter is used to convert the dust opacity obtained from the fit in a total mass. Since the empirical value comes from observations in the galactic plane, it is valid only in environments relatively similar to the ISM, specially in terms of dust-to-gass mass ratio which typical value is 0.01. Our models with the Meudon PDR code also take into account this standard value of the dust-to-gas mass ratio and the comparison with the observations have shown that this value is most likely correct (see Appendix B).

Another estimation of the molecular mass, independent of the dust-to-gas mass ratio can be obtained using emission of molecular isotopologues. Sahai et al. (2012) have estimated, using their measurement of 12CO and 13CO (3-2) line fluxes, a mass of $0.35$ M☉. With their observations, we use a slightly different method (see details in App. D.1) and obtained a value of $0.187$ M☉ under the assumption that the 13CO line is optically thin while the 12CO is optically thick. Calculations with the online non-LTE molecular radiative transfer code RADEX (van der Tak et al. 2007) confirms that the 12CO line is optically thick but indicates that the 13CO starts to be thick considering the size and the density obtained from our models. Hence, this value of $0.187$ M☉ is a lower limit.

Overall, estimates of masses are consistent but it is possible that values based on the dust emission are slightly overestimated. Using a dust-to-gas mass ratio a bit higher than the typical value of 0.01, in agreement with what is often found for protoplanetary disks (Williams & Best 2014), could correct the small differences.

D.3 Mass-loss rate at ionisation front

Recombinations of electrons and protons in the HII region give rise to the H\textalpha emission at the ionisation front with an intensity depending on the electron density $n_{\mathrm{e}}$ (see Tielens 2005, chap. 8). From the observed emission of $4\times 10^{-7}$ W m-2 sr-1 in the HST/ACS/F658N filter (including H\textalpha and [NII] 6583 Å lines), we have derived an upper limit for the electron density of $n_{\mathrm{e}}\leq 680\leavevmode\nobreak\ \mathrm{cm^{-3}}$ when considering that all the flux is coming from the H\textalpha line, which is probably close to the reality. Assuming a constant electronic density at the ionisation front all around the envelope, a shape of an ellipsoid of revolution with a surface $S$, and that the gas is escaping at the sound speed $c_{\mathrm{S}}$ calculated here for a temperature of $10^{4}$ K, the mass-loss rate is

[TABLE]

where $\mu m_{\mathrm{H}}$ is the mean particle weight. With our values, this gives a mass-loss rate of the order of $2\times 10^{\mathrm{-6}}$ M${}_{\sun}/$year for 105-600.