Widening of Protostellar Outflows: an Infrared Outflow Survey in Low Luminosity Objects

TL;DR

This study investigates the widening of protostellar outflows in low luminosity objects to understand their evolutionary stages, using infrared observations and radiative transfer modeling to measure outflow cavity angles.

Contribution

It provides the first systematic measurement of outflow opening angles in low luminosity protostars and links these angles to evolutionary indicators.

Findings

Outflows detected in 8 of 20 LLOs at Ks band.

Outflow opening angles increase with bolometric temperature.

At least 40% of LLOs are young Class 0 objects.

Abstract

We present an outflow survey toward 20 Low Luminosity Objects (LLOs), namely protostars with an internal luminosity lower than 0.2 Lsun. Although a number of studies have reported the properties of individual LLOs, the reasons for their low luminosity remain uncertain. To answer this question, we need to know the evolutionary status of LLOs. Protostellar outflows are found to widen as their parent cores evolve, and therefore, the outflow opening angle could be used as an evolutionary indicator. The infrared scattered light escapes out through the outflow cavity and highlights the cavity wall, giving us the opportunity to measure the outflow opening angle. Using the Canada-France-Hawaii Telescope, we detected outflows toward eight LLOs out of 20 at Ks band, and based on archival Spitzer IRAC1 images, we added four outflow-driving sources from the remaining 12 sources. By fitting these images with radiative transfer models, we derive the outflow opening angles and inclination angles. To study the widening of outflow cavities, we compare our sample with the young stellar objects from Arce & Sargent 2006 and Velusamy et al. 2014 in the plot of opening angle versus bolometric temperature taken as an evolutionary indicator.Our LLO targets match well the trend of increasing opening angle with bolometric temperature reported by Arce & Sargent and are broadly consistent with that reported by Velusamy et al., suggesting that the opening angle could be a good evolutionary indicator for LLOs. Accordingly, we conclude that at least 40% of the outflow-driving LLOs in our sample are young Class 0 objects.

Tables

Table 1. Table 1 Source properties

Source	Other name	R.A.	Dec	$L_{\text{int}}$	$T_{\text{bol}}$	d	Ks-band	H2 jet	Op. angle	Inc. angle	Ref.
				(L⊙)	(K)	pc	image		(∘)	∘)
DCE 001	IRAM 04191	04:21:56.88	+15:29:46.0	0.05	28 (1)	140	U∗	Y	25${}_{-7.5}{}^{}{}_{}{}^{+7.5}$	20${}_{-17.5}{}^{}{}_{}{}^{+2.5}$	1, 2
DCE 004	L1521F	04:28:38.90	+26:51:35.6	0.03	20 (1)	140	U∗	N	60${}_{-7.5}{}^{}{}_{}{}^{+2.5}$	15${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	4, 11
DCE 024	CB130-1-IRS1	18:16:16.39	-02:32:37.7	0.07	56 (7)	270	O	N	15${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	55${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	9
DCE 025	L328-IRS	18:16:59.47	-18:02:30.5	0.07	68 (7)	270	O	N	25${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	50${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	6, 10
DCE 031	L673-7	19:21:34.82	+11:21:23.4	0.04	24 (4)	300	U	N	-	-	7
DCE 032	L1148-IRS	20:40:56.66	+67:23:04.9	0.09	145 (6)	325	E	N	-	-	8
DCE 038	L1014-IRS	21:24:07.60	+49:59:08.9	0.09	67 (14)	250	O	N	65${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	3, 5
DCE 063		03:27:38.26	+30:13:58.8	0.2	199 (26)	250	E	N	-	-
DCE 064		03:28:32.57	+31:11:05.3	0.03	65 (7)	250	O	N	55${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 065		03:28:39.10	+31:06:01.8	0.02	29 (1)	250	U	N	-	-
DCE 078		03:29:23.47	+31:33:29.5	0.2	60 (5)	250	O	N	40${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 081		03:30:32.69	+30:26:26.5	0.06	33 (2)	250	U	N	-	-
DCE 090		03:32:29.18	+31:02:40.9	0.2	114 (13)	250	O	N	45${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	25${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 092		03:33:14.38	+31:07:10.9	0.14	47 (5)	250	U∗	N	25${}_{-2.5}{}^{}{}_{}{}^{+12.5}$	15${}_{-12.5}{}^{}{}_{}{}^{+7.5}$
DCE 107		03:44:02.40	+32:02:04.9	0.15	77 (3)	250	O	Y	55${}_{-27.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 109		03:44:21.36	+31:59:32.6	0.11	348 (14)	250	E	Y	-	-
DCE 145		15:40:51.62	-34:21:04.7	0.03	DCE 145 is identified as a galaxy based on the Ks-band image
DCE 181		16:26:48.48	-24:28:38.6	0.05	429 (23)	125	P	N	-	-
DCE 182	LFAM 26	16:27:05.23	-24:36:29.5	0.15	105 (6)	125	O	N	70${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 185	IRAS 16253-2429	16:28:21.60	-24:36:23.4	0.09	31 (1)	125	U∗	Y	35${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+2.5}$

Table 2. Table 2 Results of fits to Ks-band and IRAC1 data.

	Ks-band		IRAC1		IRAC2
Source	Op. angle	Inc. angle	Op. angle	Inc. angle	Op. angle	Inc. angle
	(deg	(deg)	(deg)	(deg)	(deg)	(deg)
DCE 004	-	-	60${}_{-7.5}{}^{}{}_{}{}^{+2.5}$	15${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	70${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	10${}_{-2.5}{}^{}{}_{}{}^{+2.5}$
DCE 064	55${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	45${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	-	-
DCE 078	40${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	40${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	20${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	-	-
DCE 090	45${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	25${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	30${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	55${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	-	-
DCE 092	-	-	25${}_{-2.5}{}^{}{}_{}{}^{+12.5}$	15${}_{-12.5}{}^{}{}_{}{}^{+7.5}$	45${}_{-17.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+7.5}$
DCE 107	55${}_{-27.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	40${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	5${}_{-2.5}{}^{}{}_{}{}^{+2.5}$	-	-

Equations

$$\log (\frac{{\theta }_{\mathrm{open}}}{\mathrm{deg}})=(0.47\pm 0.20)+(0.60\pm 0.10)\log (\frac{T_{\mathrm{bol}}}{\mathrm{yr}}).$$

$$\log (\frac{{\theta }_{\mathrm{open}}}{\mathrm{deg}})$$

$$=1.92+0.05\log (\frac{T_{\mathrm{bol}}}{\mathrm{yr}}),\mathrm{for}{T}_{\mathrm{bol}}>100K.$$

$$(\frac{{\theta }_{\mathrm{open}}}{\mathrm{deg}})=10^{-1.42\pm 0.39}\times (\frac{T_{\mathrm{bol}}}{K}{)}^{1.64\pm 0.21}.$$

Full text

1. INTRODUCTION

Protostellar outflows are commonly seen toward Young Stellar Objects (YSOs), especially at the embedded stage. Studying the outflow properties allows us to probe the nature of the central protostars indirectly. In a protostellar core, the outflow can carve out a biconical cavity that widens as the core evolves (Arce & Sargent, 2006; Offner et al., 2011; Velusamy et al., 2014, hereafter AS06 and VLT14 for the first and last references, respectively). The widening may initially originate from the precession of a collimated jet, but the small precession angles found in protostars cannot account for the large opening angles observed at a more evolved stage (Reipurth, 2000; Arce & Sargent, 2004). Arce & Sargent (2004), therefore, suggest that the widening of the outflow cavity is likely produced by the widening of the stellar wind from the central YSO.

The near-infrared (NIR) continuum emission traces the scattered light that escapes through the outflow cavity, and has thus offered us the opportunity to study the outflow structures over the past decades (Lucas & Roche, 1997; Padgett et al., 1999; Eisner et al., 2005; Terebey et al., 2006; Stark et al., 2006; Tobin et al., 2007; Seale & Looney, 2008; Velusamy et al., 2014). For a given outflow opening angle and a given inclination angle, the scattered light image can be modeled with radiative transfer codes (Whitney et al., 2003a, b; Robitaille et al., 2007). Therefore, comparing the synthetic images with observed NIR images allows us to derive the outflow opening angles, and thus study their evolution.

Very Low Luminosity Objects (VeLLOs), first discovered by the Spitzer Space Telescope (Young et al., 2004), are defined as embedded protostars with an internal luminosity $L_{\rm int}$ $<$ 0.1 $L_{\odot}$ (Di Francesco et al., 2007); the internal luminosity is the total luminosity of the central protostar and circumstellar disk and, at the early stage, is likely dominated by the photospheric and accretion luminosities. Using the data from the Spitzer Legacy Project “From Molecular Cores to Planet Forming Disks” (c2d; Evans et al. 2003, 2009), Dunham et al. (2008) identified 15 VeLLOs in the Perseus, Ophiuchus, Serpens, Lupus, and Chamaeleon molecular clouds plus 82 regions that contain 95 small, dense cores. During the last decade, studies of individual VeLLOs found that their properties vary much from one object to the other (IRAM 04191: André et al. 1999; Belloche et al. 2002; Dunham et al. 2006, L1014: Bourke et al. 2005; Huard et al. 2006, L1521F: Bourke et al. 2006; Takahashi et al. 2013, Cha-MMS1: Belloche et al. 2006, Tsitali et al. 2013, L328: Lee et al. 2009, 2013, L673-7: Dunham et al. 2010a, CB130: Kim et al. 2011, L1148: Kauffmann et al. 2011, IC 348-SMM2E: Palau et al. 2014, IRAS 16253: Hsieh et al. 2016). Based on these studies, Dunham et al. (2014) summarize three interpretations of the low luminosity, suggesting that VeLLOs can be (1) very low mass protostars, (2) extremely young protostars, or (3) protostars in a quiescent phase of the episodic accretion process, in which a protostar is at a quiescent accretion phase for most of the time and accretion bursts occasionally occur to deliver material onto the central protostar (Kenyon & Hartmann, 1995; Lee, 2007; Dunham et al., 2010b; Dunham & Vorobyov, 2012; Jørgensen et al., 2015; Kim et al., 2016). Schwarz et al. (2012) conducted an outflow survey of VeLLOs, and from the outflow forces, derived the time-averaged accretion luminosities; they found that these time-averaged accretion luminosities are higher than their current internal luminosities, suggesting that VeLLOs are in a quiescent phase of an episodic accretion process. On the other hand, using the N2D+/N2H+ abundance ratio as a chemical evolutionary indicator and the line width as a dynamical evolutionary indicator, Hsieh et al. (2015) suggest that VeLLOs tend to be young Class 0 protostars.

In this paper, we aim at studying the evolution of the outflow opening angle in protostellar objects at early stage. We describe the sample and the observations in Section 2. In Section 3, we report the observational results and describe how we derive the outflow opening angles using the radiative models from Whitney et al. (2003a, b). In Section 4, we discuss the correlation of the derived outflow opening angle with bolometric temperature ($T_{\rm bol}$, in comparison with AS06 and VLT14) and with the ratio of bolometric to submillimeter luminosity ($L_{\rm bol}/L_{\rm smm}$), both taken as evolutionary indicators. Finally, we summarize these results in Section 5.

2. OBSERVATIONS

2.1. Sample

We selected the faintest targets from the catalog in Dunham et al. (2008) which includes 50 low-luminosity protostars with $L_{\rm int}\leq 1.0\ L_{\odot}$ (15 have $L_{\rm int}\leq 0.1\ L_{\odot}$) in five nearby molecular clouds ($d<400$ pc) mapped by the c2d team (Evans et al., 2003, 2009). Hereafter we define the source name as the initials of the first three authors followed by the source number in Dunham et al. (2008), e.g., DCE 185. Our sample includes 13 VeLLOs out of 15 from Dunham et al. (2008), in which DCE 018 and 161 are excluded due to their low elevation ($<$35°) from the Canada-France-Hawaii Telescope (CFHT) on Mauna Kea, Hawaii. In addition to these 13 VeLLOs, we select 7 targets out of 9 in Dunham’s catalog with $0.1\ L_{\odot}\leqslant L_{\rm int}\leqslant 0.2\ L_{\odot}$, namely Low Luminosity Objects (hereafter LLOs) in order to enlarge our sample. Table 1 lists the internal luminosities and the bolometric temperatures ($T_{\rm bol}$) that separate the targets into 13 Class 0 objects ($T_{\rm bol}<$ 70 K) and 7 Class I objects (70 K $\leqslant T_{\rm bol}\leqslant$ 650 K). We note that this sample includes all 15 objects of our previous study of the envelope properties (Hsieh et al., 2015).

2.2. Ks-band images

The observations were carried out using the Wide-field InfraRed Camera (WIRCam) at the 3.6 m CFHT in the Ks broad-band filter in March and April 2010. The WIRCam has a field of view (FOV) of $20\arcmin\times 20\arcmin$ with a sampling of 0$\farcs$3 per pixel. The Ks-band filter has a bandwidth of 0.325 $\mu$m and a central wavelength of 2.146 $\mu$m. Every target was observed with a cycled dither pattern of five positions. The exposure time per image was 25 s, and the total exposure time was 300–900 s per source depending on its distance. Saturated sources were later observed with several exposures of 5 s in July and October 2011. The seeing during our observations was about 0$\farcs$7 from 0$\farcs$5 to 1$\farcs$2. The data were processed with the CFHT WIRCam standard pipeline. The data reduction was done with the TERAPIX111http://terapix.iap.fr software, that corrects the astrometry and image distortion against the Two Micron All Sky Survey (2MASS) catalog. The image flux levels were done by comparing with the 2MASS catalog.

2.3. H2 images

The H2 observations were carried out with the WIRCam at the CFHT in July and October 2011. The H2 narrow-band filter has a bandwidth of 0.032 $\mu$m and a central wavelength of 2.122 $\mu$m. The seeing was about 0$\farcs$6 (from 0$\farcs$5 to 1$\farcs$0) during these observations. The targets were observed with five dithered exposures. An individual exposure of 200 s was taken in each frame and was repeated 11–16 times for each object. Again, in order to avoid saturation issues, we took several shorter exposures of 13 s for the bright objects at NIR wavelength. The data calibration and reduction are the same as for the Ks-band observations (see Section 2.2).

2.4. Archival Spitzer data

We use the archival Spitzer InfraRed Array Camera (IRAC) images obtained by the c2d222http://peggysue.as.utexas.edu/SIRTF/ team (Evans et al., 2003, 2009). The IRAC channel 1 and channel 2 (hereafter IRAC1 and IRAC2, respectively) have central wavelength of 3.6 $\mu$m and 4.5 $\mu$m. The spatial resolutions of IRAC1 and IRAC2 are 1$\farcs$66 and 1$\farcs$72, and the pixel sizes are 1$\farcs$22 and 1$\farcs$21, respectively. These data were processed through the c2d standard pipeline, as described in the c2d data delivery document (Evans et al., 2007).

3. RESULTS AND ANALYSIS

3.1. Observed images

3.1.1 Ks-band and IRAC images

Our Ks-band observations reveal a number of sources driving protostellar outflows. Figure 1 shows the Ks-band images of the 20 selected LLOs. Thirteen objects are detected (see Table 1). Based on the morphology in the Ks-band image, we identify DCE 145 as a background galaxy misidentified as a VeLLO and exclude it from our sample. Out of the remaining 12 detected objects, eight sources show a conical or biconical structure (see Figure 1 and Table 1). The other detected sources show extended emission except for DCE 181. Although the extended emission may come from a nearly pole-on outflow, we cannot rule out the possibility that it is an infrared nebula or a background galaxy. In addition, we cannot derive the inclination angle of a nearly pole-on outflow because the inclination angle is estimated based on the contrast between the red- and blue-shifted emissions. Thus, we exclude those extended objects in our analysis. As a result, we identify eight outflow-driving LLOs on the basis of the Ks-band data.

Because the Ks-band images have a higher spatial resolution than the IRAC1 images, they allow us to explore the outflow structures better. However, the scattered light at the short wavelength (2.146 $\mu$m) could be greatly attenuated by the circumstellar envelope for deeply embedded protostars. Since the Spitzer IRAC1 images at a longer wavelength are less affected by extinction, they reveal outflow cavities in four additional embedded objects in our sample (DCE 004/L1521F, DCE 001/IRAM 04191, DCE 185/IRAS 16253, and DCE 092). In addition to extinction, the non-detections in the Ks band can also be explained by nearly edge-on configurations, as indicated by the modeling of the infrared images (see Section 3.2.2). As a result, outflow cavities are found in 12 LLOs out of 19 on the basis of the Ks-band or IRAC1 images.

To study the outflows, we compare their morphologies at NIR wavelength with that in CO emission reported in the literature. In the plane of the sky, the outflow orientations in the NIR observations are approximately consistent with that from the CO observations in DCE 004 (L1521F: Takahashi et al. 2012), DCE 025 (L328-IRS: Lee et al. 2013), DCE 038 (L1014-IRS: Bourke et al. 2005), DCE 185 (IRAS 16253: Stanke et al. 2006; Hsieh et al. 2016), DCE 092, and DCE 078 (M. Hiramatsu, private communication). However, the NIR and CO observations could sometimes trace different components of an outflow. Toward DCE 038, the CO (2–1) map (see Figure 1 in Bourke et al. 2005) reveals a compact ($\sim$5″) bipolar outflow, but our Ks-band image shows a more extended structure ($\sim$10″) with a large opening angle which is consistent with the H-band and Ks-band images in Huard et al. (2006). Although the CO observations trace large-scale outflows ($\gtrsim$100″) in DCE 001 (IRAM 04191: André et al. 1999, Lee et al. 2002) and DCE 025 (L328-IRS: Lee et al. 2013), the infrared images reveal only the inner regions ($\lesssim$10″) of the outflows. In addition, the position angle found in the IRAC1 image of DCE 001 differs by $\sim$30 degree from the position angle in the CO observations (André et al., 1999; Lee et al., 2002). Besides, the opening angle derived from the CO observations (45$\arcdeg$) is much larger than the angle derived from our fitting (25$\arcdeg$) (see Section 3.2.2). The reason for this discrepancy remains unclear, and high-angular-resolution, high-sensitivity data at NIR wavelength would be required to understand it.

Through NIR scattered light, the outflow cavities show quite different morphologies. The cavity size varies from few hundreds to few thousands au. Most of the bipolar cavities are asymmetric, being detected (or stronger) in one side only as commonly seen at NIR wavelength due to inclination effects; the undetected/weak lobes are believed to be the red-shifted lobes hidden behind the foreground envelopes. As a matter of fact, the asymmetries in the infrared images are consistent with the blue/red-shifted lobes detected by CO observations in DCE 001 (André et al., 1999) and DCE 038 (Bourke et al., 2005). CO blue-shifted and red-shifted emissions appear in both sides of the outflow in DCE 004 (Takahashi et al., 2013) and DCE 025 (Lee et al., 2013), suggesting a nearly edge-on inclination and a probably large opening angle; this result is compatible with the NIR observations. However, for DCE 185, the brighter NIR lobe in the south-west overlaps with the red-shifted lobe found by multitransition CO observations (Hsieh et al., 2016). We speculate that this inconsistency is due to an intrinsic property; instead of inclination effects, the south-west lobe may contain more scattered light because of heating by the UV-radiation from the stronger south-west H2 jet (Figure 2).

This discussion reveals that the NIR asymmetric bipolar cavities provide clues to estimate the inclination angle of the outflows. In most cases, the bipolar cavities would appear more symmetric in the systems with smaller inclination angles (with respect to the plane of the sky) and would be more asymmetric in those with larger inclination angles. Thus, we use this property to derive the inclination angle by modeling the observed images (see Section 3.2)

3.1.2 H2 images

Figure 3 shows the H2 images toward the 19 LLOs (DCE 145 has been excluded, see Section 3.1.1). The H2 observations were performed with a narrow filter but with a long integration time, allowing them to trace the continuum emission as well as the Ks-band image. The H2 images resembling the Ks-band ones are likely dominated by the continuum emission (Figure 3). Hereafter, we call Figure 3 as “H2 image(s)” when it actually consist of H2 emission and continuum emission. We note that the Ks broad-band filter includes the H2 line, but with a 10 times larger filter width corresponding to $\sim$45000 km s*-1*. Thus, the H2 emission is greatly diluted in the Ks-band image.

Comparing with the Ks-band images, our H2 survey reveals a robust jet only toward DCE 185 (IRAS 16253, Figure 2) as in the 2.12 $\mu$m H2 image of Khanzadyan et al. (2004) and in the mid-infrared InfraRed Spectrometer (IRS) H2 image of Barsony et al. (2010). The H2 jet of DCE 185 shows a prominent “S-shaped” symmetry around the central source (Figure 2). It is believed to originate from the tidal interactions between the disk where the jet originates and a companion in a noncoplanar orbit (Hsieh et al., 2016).

Besides, we identify three sources with marginal H2-jet detections by eye. Because the H2 image can be dominated by the continuum emission, we search for jet-like structures that are seen in the H2 image (Figure 3) but not in the Ks-band image (Figure 1). As a result, we find three H2-jet driving candidates: DCE 001, 109, and 107. These sources show very different morphologies of jets: (1) DCE 001 shows very weak H2 emission roughly along the outflow direction seen in the CO (2–1) map from the IRAM 30 m telescope (André et al., 1999). (2) In DCE 109, a collimated H2 emission appears from the central star to the north and may extend to a clear H2 jet knot. (3) We found two H2 knots in the north-west direction that are likely driven by DCE 107 along a collimated jet. Because DCE 001, 109, and 107 have only marginal detections of H2 jets, we call these three objects “jet-driving candidates”.

Statistically, our H2 jet survey has a low detection rate: only one prominent H2 jet-driving source and three jet-driving candidates out of 19 LLOs. These three candidates, if real at all, show very weak H2 emission. Because the H2 line emission usually traces high-velocity gas ($>$100 km s*-1*) that produces high-energy photons exciting H2 (Wolfire & Königl, 1991), this result implies that most of the LLOs drive low-velocity outflows rather than high-velocity jets. Therefore, the low detection rate of H2 jets is consistent with a scenario in which LLOs are at a quiescent accretion phase and thus drive weak outflows.

3.2. Infrared image modeling

Measuring the outflow opening angle requires a good knowledge of the inclination angle because, due to inclination effects, the apparent opening angle seen in the image is in fact larger than the true opening angle. Inclination angles can be estimated based on the intensity ratio between the blue-shifted and red-shifted lobes. Considering a spherically-symmetric envelope, the scattered light from the red-shifted lobe is more attenuated, by the thick foreground envelope, than that from the blue-shifted lobe. Therefore, a large inclination angle yields a large difference in brightness between the blue-shifted and red-shifted outflows, which allows us to derive the inclination angle, and in turn the opening angle.

We fit the NIR images (Ks-band or IRAC1) with the radiative transfer model of Whitney et al. (2003a) for the 12 sources with outflow detections (see Section 3.1.1 and Table 1). Note that although we tried to fit the three extended sources (Table 1), the best-fit results were very different from the observations (see Section 3.1.1). Thus, we exclude these three sources in our analysis. We remind the reader that the IRAC1 images have a higher sensitivity but a much lower angular resolution (see Section 3.1.1). Thus, the IRAC1 images enable us to probe outflows from more embedded objects, and the Ks-band images provide better constraints on the opening angles.

3.2.1 Setup of the models

Whitney et al. (2003a) provide a radiative transfer code to model the scattered light images of an outflow cavity (hereafter Whitney’s code/model). This code is based on the Monte Carlo radiative equilibrium routine developed by Bjorkman & Wood (2001), which calculates the radiative transfer in a three-dimensional spherical-polar grid. Whitney’s code models a spherically-symmetric envelope with biconical cavities carved out by a bipolar outflow viewed at different inclinations. The default cavity shape follows $z\propto\sqrt{x^{2}+y^{2}}^{b}$ where $x$, $y$, and $z$ are the Cartesian coordinates and $b$ is the power-law exponent of 1.5. Using this code, Whitney et al. (2003b) present an evolutionary sequence of models for a low mass protostar at stages of Class 0, Late 0, I, Late I, II, and III. Because our targets have bolometric temperatures ranging from 20 K to 114 K (Table 1), we take the Class Late 0 model as the template model and then vary the opening angle and inclination angle to fit the observed Ks-band and IRAC1 images.

Here, we discuss how the parameter set in Whitney’s code may affect the modeled images. Whitney’s code categorizes several tens of parameters into three groups: central star, disk, and envelope. We speculate that the central star and disk may only weakly affect the large-scale outflow cavities at hundreds to thousands au, when the disk of a Class Late 0 source (Whitney et al., 2003b) has an outer radius of 50 au ($\sim$0$\farcs$2–0$\farcs$4 in our sample). In addition, to remove the influence of the star and disk emission, we mask the central source in a circular region in our fitting process (see Section 3.2.2). Furthermore, we find that the internal luminosity (associated to parameters of star and disk) likely affects only the image brightness scale of the cavity but not the distribution or structure. Therefore, we use a linear function to scale the intensity map (see Section 3.2.2), which makes our fitting results independent of the internal luminosity. On the other hand, we test whether or not the envelope properties could change the brightness distributions of the outflow cavities. We find that the density profiles in the cavity can significantly change the intensity scale but not the distribution. We also find that the mass infall rate ($\dot{M}_{\rm env}$) can alter the intensity distribution but does not significantly affect the opening angle, and for the Class Late 0 objects, it is set to 10*-5* $M_{\odot}\rm~{}yr^{-1}$ (Whitney et al., 2003b).

3.2.2 Fitting the opening angles and inclination angles

There are several steps in our fitting process. Since the luminosity of the central source can change the image brightness, we linearly scale the intensity of the modeled image to fit the observations. To remove the effects of field stars and the central object (star$+$disk), we mask these sources with adequate circles determined by eye. Using Whitney’s code, we construct a grid of models that vary opening angle ($\theta_{\rm open}$), inclination angle ($\theta_{\rm inc}$), and position angle (P.A.) with a cell size of 5$\arcdeg$ for all dimensions. Then, we derive the $\chi^{2}$ values between the observed image and each model, and as a result, we obtain a $\chi^{2}$ distribution in the three-dimensional grid. Then, we find the best-fit opening angle and inclination angle, and taking a 99.7% confidence level with five free parameters ($\theta_{\rm open}$, $\theta_{\rm inc}$, P.A., and two in the linear function) ($\Delta\chi^{2}>18.2$, Press et al. 1992), we obtain its upper and lower limits (Table 1). Figure 4 shows the images of the best-fit models and the corresponding observed images.

We discuss here two special cases, DCE 185 and 092, with arbitrary decisions in our fitting process. For DCE 185, as the south-west lobe is brighter than the north-east lobe in the IRAC1 image (Figure 2), this brightness asymmetry is in conflict with the blue- and red-shifted outflows identified by CO observations (Stark et al., 2006; van der Marel et al., 2013; Hsieh et al., 2016). Thus, we speculate that this asymmetric bipolar cavity is due to an intrinsic property rather than inclination effects (see Section 3.1.1). In the fitting, we restrict the red-shifted lobe to be in the south-west as seen in the CO observations, and since the model assumes the asymmetry comes from inclination effects, we obtain an inclination angle of 5°, the lower limit of our model grid. This underestimate of the inclination angle yields an overestimate of the opening angle. Therefore, we take the derived opening angle of DCE 185 as an upper limit for later analysis. Besides, the $\chi^{2}$ distribution of DCE 185 shows two local minima at opening angles of 125$\arcdeg$ and 35$\arcdeg$. We eliminate the larger opening angle, since it is most likely caused by contamination of diffuse cloud emission. For DCE 092, the $\chi^{2}$ distribution also has two local minima located at opening angles of 100$\arcdeg$ and 35$\arcdeg$. We eliminate the larger opening angle because it is likely affected by the nearby bright sources.

To check our fitting results, we fit the IRAC1 images in four sources (DCE 064, 107, 090, and 078) that have outflow detections at IRAC1 in addition to the Ks-band, and the IRAC2 images of the following two sources: (1) DCE 004, which is an outlier in the $T_{\rm bol}-\theta_{\rm open}$ plot (see Section 4.1 and Figure 6), and (2) DCE 092, for which the inclination angle derived from the fitting of its IRAC1 image has a large uncertainty. Comparing the fitting results at two wavelengths, we obtain consistent inclination angles and similar opening angles (difference $\leqslant$ 10° or within the errors), except for DCE 090 (Table 2 and Figure 5). In DCE 090, we suspect that the discrepancy between $\theta_{\rm inc,IRAC1}=55\arcdeg$ and $\theta_{\rm inc,Ks}=25\arcdeg$ comes from an unreliable fitting in IRAC1 due to the contamination from the bright central object. Because the Ks-band image has a better resolution, we believe that the result from the Ks-band is more reliable for DCE 090. Overall, we conclude that fitting the images at different wavelengths gives comparable results. As a result, we list the best-fit results from the relatively shorter wavelengths in Table 1 and use these values in the following analysis.

3.2.3 Fitting with streamline cavity

In addition to the curved cavity mentioned above, Whitney’s code provides an alternative cavity shape namely streamline which is conical on large scales and might be carved out by precessing jets (Whitney et al., 2003a). To study the origin of the outflow widening, we also fit the observed images with the model of streamline cavity and compare the results with those obtained assuming a curved cavity. We find that, except for DCE 001 and DCE 024 that are better fitted with the streamline cavity, all other objects are better fitted with the curved cavity above a confidence level of 99.7%. DCE 001 and DCE 024 have the smallest opening angles and relatively low bolometric temperatures among our targets, suggesting that they are younger than other sources. Although, based on this, one could be tempted to conclude that the outflow cavity is carved by precessing jets at the earliest stage, a sample of only two sources is not large enough to be statistically meaningful. In addition, the large-scale CO outflows in DCE 001 do not show the streamline shape (André et al., 1999). Therefore, we lack evidence to support the fact that the outflow widening is caused by carving of precessing jets.

4. DISCUSSION

4.1. Outflow opening angle versus bolometric temperature

Here, we compare our LLOs with YSOs from AS06 and VLT14 in the plot of bolometric temperature versus opening angle (Figure 6a). To study the evolution of protostellar outflows, previous works have studied the correlation between outflow opening angle and bolometric temperature taken as an evolutionary indicator. Based on their survey of CO outflows toward YSOs, AS06 found that the outflow opening angle widens as the core evolves from Class 0 to Class III. They found a correlation between the outflow opening angle ($\theta_{\rm open}$) and the bolometric temperature ($T_{\rm bol}$) as

[TABLE]

Lately, using Spitzer IRAC images, VLT14 measured the outflow opening angles toward 31 YSOs and obtained a $T_{\rm bol}-\theta_{\rm open}$ correlation of

[TABLE]

Before we compare our sample with AS06’s and VLT14’s results, two caveats should be made. First, because both AS06 and VLT14 measured the opening angles without inclination correction, their measurements should be viewed as upper limits in comparison with ours. Second, the NIR scattered light could highlight a broader outflow cavity than CO (Ohashi et al., 1997; Tafalla & Myers, 1997; Lee et al., 2006; Launhardt et al., 2008; Arce et al., 2013), because CO may not be able to trace the full extent of the NIR reflection nebula due to its high opacity especially at the velocities close to the ambient cloud velocity (Arce et al., 2013).

We find that most LLOs fit well the $T_{\rm bol}-\theta_{\rm open}$ trend found by AS06, except for DCE 004 (L1521F) and DCE 024 (CB-130-1-IRS1). Two reasons may explain the deviation of DCE 004: (1) DCE 004 may host a binary system driving two outflows, as SMA CO (2–1) observations suggest (Takahashi et al., 2013). The possible two outflows have different axes as seen in many other cases (Offner et al., 2016; Lee et al., 2016), which can affect our opening angle estimate with NIR observation. Therefore, we consider the measured opening angle as an upper limit. (2) Our fitting result implies a nearly edge-on configuration, which may cause the low bolometric temperature (see Section 4.2). Although DCE 024 deviates from the power-law relation of AS06, it fits well the new power-law index that better describes the sources with $T_{\rm bol}\lesssim 120$ K (see below).

We use a single power-law to fit our LLO sample (Figure 6a). The best fit is

[TABLE]

We have excluded DCE 004 and DCE 185 from the fitting since their fitted opening angles are considered as upper limits. The best-fit power-law index, 1.64, is much larger than the index derived by AS06. If we include the data points with $T_{\rm bol}\leqslant 120$ K from AS06, we obtain $(\frac{\theta_{\rm open}}{\rm deg})=10^{-1.02\pm 0.26}\times(\frac{T_{\rm bol}}{K})^{1.46\pm 0.14}$, which is consistent with Equation 3 within $1\sigma$. This suggests that we need a broken power law or a more complicated model to describe how the outflow cavities widen. In addition, based on this large index, we suggest that the outflow opening angle may better discriminate the evolutionary state of protostars than the bolometric temperature in a range of $T_{\rm bol}\leqslant 120~{}K$.

VLT14 used two power-law components to fit their data and found a break at $T_{\rm bol}\approxeq 100~{}K$. However, the steeper power law at the lower $T_{\rm bol}$ end (Equation 2a) is quite different from our fitting result (Equation 3). We speculate that this discrepancy partially comes from inclination effects that were not taken into account by VLT14. To correct for this, we assume that all sources of VLT14 (as well as AS06) have an inclination angle of 32°7 (with respect to the plane of the sky) and we estimate their inclination-corrected opening angles (Figure 6b); the angle of 32°7 corresponds to a mean inclination angle assuming all orientations are equally favorable (Bontemps et al., 1996; Dunham et al., 2014, note that the angle 57°3 in the references is with respect to the line of sight). After the correction, the VLT14 sample is broadly consistent with our LLOs. We find a new-best fit of $\log(\frac{\theta_{\rm open}}{\rm deg})=-0.38+1.25\log(\frac{T_{\rm bol}}{\rm yr})$ for $T_{\rm bol}<100K$ for the inclination-corrected sample of VLT14, which is consistent with our LLOs (Equation 3) within $\sim$3-4$\sigma$.

We now compare the opening angles of our LLO sample with the Class 0 samples studied by AS06, Seale & Looney (2008), and VLT14. The AS06 sample includes 11 Class 0 objects, Seale & Looney’s sample has 21 Class 0 objects, and the VLT14 sample contains 20 Class 0 objects. Out of the 12 LLOs with NIR outflow detections in our sample, five ($\sim$40%) have opening angles smaller than $25\arcdeg$. Such small opening angles are rare in the Class 0 samples of AS06, Seale & Looney (2008), and VLT14 (3/11, 2/21, and 0/20, respectively). Note that the ratios remain the same after the inclination correction. Assuming that the opening angle is an evolutionary indicator, we conclude that our sample of outflow-driving LLOs contains a higher fraction of young objects than the samples of AS06, Seale & Looney (2008), and VLT14. In turn, this suggests that at least the five LLOs with $\theta_{\rm open}\leqslant 25\arcdeg$ may be very young Class 0 objects.

4.2. Approximate correction for inclination effects on $T_{\rm bol}$

Figure 6a shows that the LLOs below the best-fit line have relatively large inclination angles while the sources above have relatively small inclination angles. Because the inclination angle can affect the Spectral Energy Distribution (SED), and in turn the bolometric temperature (Whitney et al., 2003a, b; Robitaille et al., 2006; Crapsi et al., 2008), this geometrical effect may cause the LLOs to deviate from the best-fit relation.

To investigate whether the inclination effects cause the dispersion in Figure 6a, we correct the bolometric temperatures for inclination effects using the SED fitting tool of Robitaille et al. (2006, 2007). This SED fitting tool, constructed with the radiative transfer code of Whitney et al. (2003a), consists of 20,000 YSO models with SEDs computed at 10 viewing angles for each model, resulting in 200,000 SEDs in total. Since the SEDs’ uncertainty grows rapidly at wavelengths $\gtrsim$100 $\mu$m and becomes very high at wavelengths $\gtrsim$600 $\mu$m (see Figure 1 in Robitaille et al., 2006), we ignore the SED data points at wavelengths $>$850 $\mu$m, and for the sources lacking 850 $\mu$m data, interpolate a flux at 850 $\mu$m logarithmically. From the processed SEDs, we calculate the bolometric temperatures (Table 1). The derived temperatures are consistent with the results of Dunham et al. (2008) within the uncertainties except for DCE 031 (DCE 031 has no outflow detection such that it will not affect the following analysis.). We fit the processed SEDs with the YSO models of Robitaille et al. (2006, 2007). Then, we identify the best-fit model as the one with the lowest $\chi^{2}$ and with the consistent $\theta_{\rm open}$ and $\theta_{\rm inc}$ from our image fits (see Section 3.2.2) within a difference $<$10°. To get rid of inclination effects on $T_{\rm bol}$, we define an inclination-independent bolometric temperature as the bolometric temperature of the synthetic SED of the best-fit model computed for an arbitrary inclination angle of 27°. This inclination angle corresponds to about the averaged value in our sample (Table 1). We remind the readers that our survey at NIR wavelength cannot identify outflows at a nearly pole-on configuration. Therefore, a small inclination angle is reasonable. We calculate the bolometric temperatures from the best-fit SED and from the synthetic SED, namely best-fit $T_{\rm bol}$ and synthetic $T_{\rm bol}$, respectively. With the same inclination, we can use the synthetic $T_{\rm bol}$ as an evolutionary indicator without inclination effects. Figures 7a and 7b show the $T_{\rm bol}-\theta_{\rm open}$ relations with the best-fit $T_{\rm bol}$ and synthetic $T_{\rm bol}$, respectively, in which, for comparison, we assign the same uncertainty as the observed $T_{\rm bol}$ in percentage. By fitting a single power-law as Equation (3), we find that the $T_{\rm bol}-\theta_{\rm open}$ correlation with the synthetic $T_{\rm bol}$ has a smaller $\chi^{2}$ (see Figure 7). This result supports our speculation that the inclination effects may cause the deviation from the best-fit $T_{\rm bol}-\theta_{\rm open}$ relation.

There are three caveats about the inclination correction. First, although we use the image fits to constrain $\theta_{\rm open}$ and $\theta_{\rm inc}$, the SED fit may still be degenerated, preventing us from determining the physical parameters. Thus, our inclination correction on $T_{\rm bol}$ should be considered as a rough approximation rather than a precise determination. Second, the SED models of Robitaille et al. (2006) do not include brown dwarfs while they are used to interpret the low luminosity of VeLLOs in the literature (Bourke et al., 2005; Kauffmann et al., 2011; Lee et al., 2013; Palau et al., 2014). Third, the synthetic $T_{\rm bol}-\theta_{\rm open}$ correlation is probably a self-consistent result produced artificially by the radiative transfer code (Whitney et al., 2003a; Robitaille et al., 2006, 2007) because both the synthetic $T_{\rm bol}$ and $\theta_{\rm open}$ correspond to the model rather than the observation. Despite these caveats, the synthetic $T_{\rm bol}-\theta_{\rm open}$ and the observed $T_{\rm bol}-\theta_{\rm open}$ have very similar power-law fits (see equations in Figures 6 and 7b). This implies that the radiative transfer code reproduces well the observed $T_{\rm bol}-\theta_{\rm open}$ correlation at least for these LLOs. This discussion tentatively suggests that the tighter synthetic $T_{\rm bol}-\theta_{\rm open}$ correlation is reliable.

To further test the inclination correction on $T_{\rm bol}$, we compare the best-fit $T_{\rm bol}$ and the synthetic $T_{\rm bol}$ with the ratio of bolometric to submillimeter luminosity ($L_{\rm bol}/L_{\rm smm}$) taken as a better evolutionary indicator. We exclude DCE 038 because of its problematic measurement of $L_{\rm bol}/L_{\rm smm}$ (see Section 4.3). We find that $L_{\rm bol}/L_{\rm smm}$ is slightly better correlated to the synthetic $T_{\rm bol}$ with $r=0.57$ and $\rho=0.69$ (p-values of 0.11 and 0.04, respectively) than to the best-fit $T_{\rm bol}$ with $r=0.09$ and $\rho=0.50$ (p-values of 0.81 and 0.17, respectively). Although the sample is small, this result suggests that the synthetic $T_{\rm bol}$ could better reflect the evolutionary stage than the best-fit $T_{\rm bol}$.

4.3. Outflow opening angles versus $L_{\rm bol}/L_{\rm smm}$ and other parameters from SED

To study what physical conditions may affect the outflow opening angle, we compare the latter to four parameters from Dunham et al. (2008): the internal luminosity ($L_{\rm int}$), the bolometric luminosity ($L_{\rm bol}$), the ratio of bolometric to submillimeter luminosity ($L_{\rm bol}/L_{\rm smm}$), and the submillimeter luminosity ($L_{\rm smm}$). We exclude DCE 185 and 004 because their opening angles are considered as upper limits (see Sections 3.2.2 and 4.1). We evaluate the significance of correlation between the outflow opening angle and these parameters using Pearson’s $r$ correlation test and Spearman’s $\rho$ rank correlation test (Conover et al., 1999). Based on the correlation tests, we find no significant correlations between $\theta_{\rm open}$ and these parameters, except for $L_{\rm bol}/L_{\rm smm}$. We find a correlation coefficient $r=0.55$ for $\theta_{\rm open}-L_{\rm bol}/L_{\rm smm}$, suggesting a probability (p-value) of 10% for it being an uncorrelated/random distribution. Since $L_{\rm bol}/L_{\rm smm}$ is considered as an evolutionary indicator (André et al., 1993; Young & Evans, 2005), we conclude that the opening angle is most likely sensitive to the evolutionary status but not to other physical conditions such as envelope mass traced by $L_{\rm smm}$.

Young & Evans (2005) suggest that $L_{\rm bol}/L_{\rm smm}$ can better reflect the evolutionary status than $T_{\rm bol}$ because, in their evolutionary model, $L_{\rm bol}/L_{\rm smm}$ is sensitive to the fraction of mass accreted onto the central object but is less affected by the initial core mass. Figure 8 shows our LLO sample in a $L_{\rm bol}/L_{\rm smm}$ versus $T_{\rm bol}$ diagram, along with the boundaries traditionally used to separate the Class 0 and I phases, either in bolometric temperature ($T_{\rm bol}=70K$ Chen et al., 1995) or in bolometric-to-submm luminosity ratio ($L_{\rm bol}/L_{\rm smm}$=200 André et al., 1993). The classification of our LLO sample into Class 0 and I objects depends somewhat on the chosen axis. Three objects out of 19 fall into the Class 0 category according to their $L_{\rm bol}/L_{\rm smm}$ ratio but into the Class I category according to their bolometric temperature. The classification of the rest of the sample is the same with both criteria. Figure 8 shows a correlation between $T_{\rm bol}$ and $L_{\rm bol}/L_{\rm smm}$ toward the 19 LLOs ($r=0.74$ and $\rho=0.72$ with both p-values $=0$), though the correlation is less significant when considering only the ten targets used to derive the $T_{\rm bol}-\theta_{\rm open}$ relation. DCE 038 seems to be an outlier in both Figures 8 and 9 with an extremely low $L_{\rm bol}/L_{\rm smm}$ ratio (3 $\pm$ 1). This source has a puzzling jump in SED with 3.2 $\pm$ 0.5 Jy at 350 $\mu$m and 21.5 $\pm$ 16.1 Jy at 450 $\mu$m (1.8 $\pm$ 0.4 Jy at 850 $\mu$m), which results in the low $L_{\rm bol}/L_{\rm smm}$ (Dunham et al., 2008) when the submillimeter luminosity is defined as the luminosity at $\lambda>350~{}\mu$m (André et al., 1993). $L_{\rm bol}/L_{\rm smm}$ would become 10 $\pm$ 1 if we ignore the flux density at 450 $\mu$m or 6 $\pm$ 3 if we replace the observed flux density at 350 $\mu$m by a linear interpolation between 70 and 450 $\mu$m. This suggests that the observed $L_{\rm bol}/L_{\rm smm}$ is underestimated. Therefore, we remove DCE 038 and, in turn, find better correlation coefficients of $r=0.74$ and $\rho=0.80$ between $T_{\rm bol}$ and $L_{\rm bol}/L_{\rm smm}$ (Figure 8). Furthermore, comparing to the outflow opening angle, we obtain correlation coefficients of $r=0.71$ and $\rho=0.77$ (p-values of 0.03 and 0.02, respectively) for $L_{\rm bol}/L_{\rm smm}-\theta_{\rm open}$ (Figure 9) and $r=0.67$ and $\rho=0.70$ (p-values of 0.05 and 0.04, respectively) for $T_{\rm bol}-\theta_{\rm open}$. Therefore, we suggest that the distribution of opening angles is better correlated with evolutionary status using $L_{\rm bol}/L_{\rm smm}$ as an age indicator than using $T_{\rm bol}$. This result implies that the outflow opening angle likely reflects the evolutionary status well and could be considered as a good evolutionary indicator.

5. SUMMARY

We conducted an outflow/jet survey of 20 Low Luminosity Objects (LLOs) using CFHT in Ks-band continuum emission and H2 line emission. From the Ks-band observations, we identify a background galaxy and eliminate it from our sample. We detect outflow cavities in eight sources out of 19 LLOs on the basis of the Ks-band continuum observations. Among the remaining 11 LLOs, the archival Spitzer IRAC1 data reveals four outflow-driving sources. We derive the outflow opening angles and the inclination angles of these 12 LLOs by fitting the observed images with the radiative transfer models of Whitney et al. (2003a). The H2 observations reveal only one LLO associated with a prominent H2 jet. Three other LLOs out of the remaining 18 have marginal detections of an H2 jet. Our main results are the following:

As we detect only one to four H2 jets out of 19 LLOs, we suggest that most LLOs do not generate strong jets or outflows. This indirectly supports the idea that these LLOs are likely at a quiescent accretion phase. 2. 2.

Our LLOs follow a trend similar to the one found by AS06 in the plot of bolometric temperature versus opening angle ($T_{\rm bol}-\theta_{\rm open}$), and are broadly consistent with the correlation reported by VLT14, after correction for the inclination. Instead of a single power-law describing the distribution, we find a larger index for the sources at an early evolutionary stage ($T_{\rm bol}$ $\lesssim$ 120 K), in agreement with VLT14. We conclude that the outflow opening angle may better trace the evolutionary stage than $T_{\rm bol}$ for $T_{\rm bol}<120$ K. 3. 3.

The LLOs located above the best-fit power-law index have relatively small inclination angles and those located below have large inclination angles. This suggests that the dispersion in the $T_{\rm bol}-\theta_{\rm open}$ plot may be in part due to inclination effects on $T_{\rm bol}$. 4. 4.

Using the outflow opening angle as an evolutionary indicator, we suggest that at least 40% of the outflow-driving LLOs in our sample are young Class 0 objects. 5. 5.

Out of the 12 targets with infrared outflow detections, ten are better fitted with a curved cavity than a streamline cavity, suggesting that the outflow cavities are not carved by precessing jets.

The authors thank Dr. Masaaki Hiramatsu for fruitful discussions and a comparison to his SMA outflow survey of LLOs in Perseus. We are grateful to Dr. Chi-Hung Yan for his assistance with the CFHT data reduction. We would like to thank Prof. Héctor Arce, Dr. Ian Stephens, and Dr. Thomas Robitaille for providing valuable discussions. The authors acknowledge the staff at CFHT for assistance with operations. The authors thank the referee for the insightful comments that improved this paper. T.H.H. and S.P.L. are thankful for the support of the Ministry of Science and Technology (MoST) of Taiwan through Grants NSC 98-2112- M-007-007-MY3, NSC 101-2119-M-007-004, MoST 102-2119-M-007-004-MY3 and MoST 105-2119-M-007-024.