THE ENIGMATIC CORE L1451-mm: A FIRST HYDROSTATIC CORE? OR A HIDDEN VeLLO?^*

Dust continuum observations at 1.2 mm were taken using MAMBO at IRAM 30 m telescope, under good weather (τ_1.2 mm = 0.1–0.2). The data reduction was carried out using MOPSI, with parameters optimized for extended sources. The observations are convolved with a 15'' Gaussian kernel, while the flux unit is in Jy per 11'' beam. The rms noise level is 1 mJy per 11'' beam, and the map for the core studied is shown in Figure 2.

Zoom In Zoom Out Reset image size

Observations were carried out with the Very Large Array (VLA) of the National Radio Astronomy Observatory on 2006 January 10 (project AA300). The NH₃ (J, K) = (1, 1) and (2, 2) inversion transitions were observed simultaneously (see Table 2 for a summary of the correlator configuration used). At this frequency the primary beam of the antennas is about 19. The array was in the compact (D) configuration, the bandwidth was 1.56 mHz, and the channel separation was 12.2 kHz (corresponding to 0.154 km s⁻¹). This configuration is centered at the main hyperfine component and it also covers the inner pair of satellite lines for NH₃ (1,1).

The bandpass and absolute flux calibrator was the quasar 0319+415 (3C84) with a calculated flux density of 10.6 Jy at 1 cm, and the phase and amplitude calibrator was 0336+323. The raw data were reduced using CASA image processing software. The signal from each baseline was inspected, and baselines showing spurious data were removed prior to imaging. The images were created using multi-scale clean (scales [8,24,72] arcsec and small-scale bias = 0.8) with a robust parameter of 0.5 and tapering the image with a 8'' Gaussian to increase the signal to noise. Each channel was cleaned separately according to the spatial distribution of the emission, using a circular beam of 8''. Table 3 lists relevant information on the maps used.

Table 3. Parameters of Interferometric Maps

Map	Array	Beama	rms
NH3 (1,1)	VLA	8'' × 8'' (0°)	3 mJy beam−1  channel−1
NH3 (2,2)	VLA	8'' × 8'' (0°)	3 mJy beam−1  channel−1
NH2D (111–101)	CARMA	92 × 76 (+723)	90 mJy beam−1  channel−1
N2H+ (1–0)	CARMA	52 × 43 (−73°)	80 mJy beam−1  channel−1
3 mm continuum	CARMA	54 × 48 (−77°)	0.5 mJy beam−1
12CO (2–1)	SMA	135 × 096 (+809)	40 mJy beam−1  channel−1
1.3 mm continuum	SMA	123 × 088 (+856)	0.5 mJy beam−1

Note. ^aSize and position angle. Position angle is measured counterclockwise from north.

Download table as:  ASCIITypeset image

Continuum observations in the 3 mm window were obtained with CARMA, a 15-element interferometer consisting of nine 6.1 m antennas and six 10.4 m antennas, between 2008 April and September. The CARMA correlator records signals in three separate bands, each with an upper and lower sideband. We configured one band for maximum bandwidth (468 MHz with 15 channels) to observe continuum emission, providing a total continuum bandwidth of 936 MHz. The remaining two bands were configured for maximum spectral resolution (1.92 MHz per band) to observe NH₂D (1₁₁–1₀₁) and N₂H⁺ (1–0) (see Table 4 for the correlator configuration summary). The six main hyperfine components of NH₂D fit in the two narrow spectral bands and six of the seven hyperfine components of N₂H⁺ (1–0) were observed, with the highest frequency (isolated) component falling outside the observed frequency range.

The field of view (half-power beam width) of the 10.4 m antennas is 66'' at the observed frequencies. Seven point mosaics were made around the center of L1451-mm in CARMA's D- and E-array configurations, giving baselines that range from 8 m to 150 m. Observations of N₂H⁺ (but not NH₂D) were also made in CARMA's C-array configuration, with projected baselines of 30–350 m. The synthesized beam sizes and position angles (measured counterclockwise from north) are 54 × 48 and −77° (continuum), 52 × 43 and −73° (N₂H⁺), 92 × 76 and 723 (NH₂D). The largest angular size to which these observations were sensitive is ∼40''.

The observing sequence for the CARMA observations was to integrate on a primary and secondary phase calibrator (3C 111 and 0336+323) for 3 minutes each and the science target for 14 minutes. In each set of observations 3C 111 was used for passband calibration and observations of Uranus were used for absolute flux calibration. Based on the repeatability of the quasar fluxes, the estimated random uncertainty in the measured fluxes is σ ≃ 5%. Radio pointing was done at the beginning of each track and pointing constants were updated at least every two hours thereafter, using either radio or optical pointing routines (Corder et al. 2010). Calibration and imaging were done using the MIRIAD data reduction package (Sault et al. 1995). Table 3 lists relevant information on the maps used.

The SMA observations were carried out at 1.3 mm (230 GHz) in both compact and extended configuration. The compact array observations were carried out on 2009 November 1, with zenith opacity at 225 GHz of ∼0.085. Quasars 3C 84 and 3C 111 were observed for gain calibration. Flux calibration was done with observations of Uranus and Ganymede. Bandpass calibration was done using observations of the quasar 3C 273. The SMA correlator covers 2 GHz bandwidth in each of the two sidebands. Each band is divided into 24 "chunks" of 104 MHz width, which can be covered by varying spectral resolution. The correlator configuration is summarized in Table 5.

The extended array observations were carried out on 2010 September 13, with zenith opacity at 225 GHz of ∼0.05. Quasars 3C 84 and 3C 111 were observed for gain calibration. Flux calibration was done with observations of Uranus and Callisto. Bandpass calibration was done using observations of the quasar 3C 454.3. The SMA correlator used the new 4 GHz bandwidth in each of the two sidebands. Each band is divided into 48 "chunks" of 104 MHz width, which can be covered by varying spectral resolution. The correlator configuration is summarized in Table 6.

Both data sets were edited and calibrated using the MIR software package¹⁰ adapted for the SMA. Imaging was performed with the MIRIAD package (Sault et al. 1995), resulting in an angular resolution of 123 × 088 P.A. = 856 (using robust weighting parameter of −2) and 135 × 096 P.A. = 809 (using robust weighting parameter of 0) for the continuum and ¹²CO (2–1), respectively. Table 3 lists relevant information on the maps used. The rms sensitivity is ≈0.5 mJy beam⁻¹ for the continuum, using both sidebands (avoiding the chunk containing the ¹²CO line), and ∼36 mJy beam⁻¹ per channel for the line ¹²CO (2–1) data. The primary beam FWHM of the SMA at these frequencies is about 55''.

The MAMBO dust continuum emission map (left panel of Figure 2) can be decomposed into a bright compact core and fainter filamentary emission. The compact core peak position is located at (α, δ) = (03:25:10.4, +30:23:56.0). The compact core within 4200 AU mass is estimated to be 0.3 M_☉, where a dust opacity per dust mass of 1.14 cm² g⁻¹ (Ossenkopf & Henning 1994), a gas-to-dust ratio of 100, and a dust temperature of 10 K are used. The MAMBO-derived mass is consistent with the mass previously estimated using Bolocam. In Figure 3, the compact core is compared to the sample of starless cores from Kauffmann et al. (2008), where the fiducial radius of 4200 AU is used to compare with previous works (e.g., Motte & André 2001) suggesting that it is more compact than most starless cores. From this comparison we can establish that L1451-mm is in fact quite compact, and therefore dense, suggesting that it is more compact than most starless cores an evolved evolutionary state (see Crapsi et al. 2005).

Zoom In Zoom Out Reset image size

A faint central source is detected in the CARMA 3 mm continuum map shown in the right panel of Figure 2. The continuum emission map is fitted by a Gaussian with a total flux of 10 mJy, while if a point source is fitted a flux of (4 ± 2) mJy is obtained. A summary of the fits to the CARMA 3 mm continuum is listed in Table 7. The CARMA continuum emission agrees with the MAMBO peak position.

Figure 4 shows the summary of the molecular line transitions observed with CARMA and VLA. In this study, we will briefly discuss the kinematics of the region and leave a more in depth study of the core in forthcoming papers (S. Schnee et al. 2011, in preparation; H. G. Arce et al. 2011, in preparation). From these observations, a centroid velocity and velocity dispersion are obtained by fitting the line profiles (see Rosolowsky et al. 2008 for details). The integrated intensity maps show the extended emission from the core where the peaks match the position of the CARMA continuum emission to within the respective beam size. The centroid velocity maps, for all three lines, show a consistent result with a clear velocity gradient. A gradient is fitted to the centroid velocity map for all three lines, with an average value of $\mathcal {G}=$ 6.1 km s⁻¹ pc⁻¹ and a −66° position angle (measured counterclockwise from north), see Table 8 for the individual fit obtained for all three maps. This velocity gradient is larger than those observed in lower angular resolution NH₃ (1,1) maps (Goodman et al. 1993) or using lower density tracers (Kirk et al. 2010), while velocity gradients of a similar magnitude are obtained with high angular resolution observations of dense gas (Curtis & Richer 2011; Tanner & Arce 2011). The velocity dispersion maps show a clear increase toward the center of the map, starting with very narrow lines (close to the thermal values) in the outer regions. The increase in velocity dispersion is more pronounced in the N₂H⁺ velocity dispersion map, and the difference can be explained by the higher angular resolution obtained in the N₂H⁺ observations.

Zoom In Zoom Out Reset image size

Table 8. Velocity Gradients Fits

Transition	$\mathcal {G}$	P.A.	v0
	(km s−1 pc−1)	(deg)	(km s−1)
N2H+	5.6 ± 0.2	−84 ± 3	3.970 ± 0.002
NH2D	8 ± 1	−83 ± 15	3.949 ± 0.009
NH3	6.24 ± 0.06	−65.8 ± 0.7	4.0042 ± 0.0006

Download table as:  ASCIITypeset image

The NH₃ (1,1) and (2,2) integrated intensity maps obtained using the VLA are shown in Figure 5, left and middle panels respectively, where all components observed are taken into account. Both lines present a peak coincident with the CARMA continuum peak and where the NH₃ (1,1) emission covers a more extended region than the NH₃ (2,2). However, since the emission NH₃ (1,1) is fairly extended, the addition of GBT data to provide the zero spacing is needed to allow a robust temperature determination. The morphology of the NH₃ and N₂H⁺ integrated intensity maps show non-flattened structures, which are drastically different from those seen in young Class 0 sources (e.g., Wiseman et al. 2001; Chiang et al. 2010; Tanner & Arce 2011). The right panel of Figure 5 presents the derived kinetic temperature obtained from the simultaneous NH₃ (1,1) and (2,2) line fit, with uncertainties in the temperature determination between 0.2 K, in the central region, up to 0.5 K in the outer regions. Surprisingly, the kinetic temperature map is quite constant, in particular, there is no evidence for an increase in temperature toward the peak continuum position.

Zoom In Zoom Out Reset image size

The visibility amplitude as a function of uv-distance for the 1.3 mm continuum is shown in Figure 6. From this figure we identify two components. An extended component that is resolved at long baselines, and a compact component that remains unresolved even at the longest baselines (indicated by the horizontal line in Figure 6). This unresolved component is commonly seen toward dense cores containing a central protostar, and it is interpreted as arising from an unresolved central disk (e.g., Jørgensen et al. 2007, 2009). However, in the case of L1451-mm there is no infrared detection of a central protostar.

Zoom In Zoom Out Reset image size

The SMA 1.3 mm continuum map is shown in Figure 7, with redshifted and blueshifted ¹²CO (2–1) emission overlaid. The dust continuum emission clearly shows the central source, also presented in Figure 6. The position of this continuum source coincides with the pixel where the CARMA 3 mm continuum peaks. The orientation of the SMA continuum emission, obtained through a fit of the visibilities and listed in Table 9, is close to the right ascension axis and clearly different from the redshifted and blueshifted ¹²CO emission. The detected ¹²CO (2–1) emission shows a large velocity dispersion, with spatially separated blueshifted and redshifted lobes (see Figure 7). Figure 8 shows the PV diagram along the gray line drawn in Figure 7, where the dashed vertical line shows the centroid velocity of the dense core, 3.94 km s⁻¹, which is consistent with the interferometric observations (see Table 8) and the NH₃ (1,1) data obtained with the GBT at 30'' (J. E. Pineda et al. 2011, in preparation).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dust mass of the compact (unresolved) emission is estimated as

$$\begin{equation} M_{1.3 \,{\rm mm}} = 1.3\, {M_{\odot }} \left(\frac{F_{1.3 \,{\rm mm}}}{1\,{\rm Jy}}\right) \left(\frac{d}{200 \,{\rm pc}}\right)^{2} \left(e^{0.36\,(30 K/T)} -1\right), \end{equation} \tag{ 4 }$$

where we assumed optically thin emission and a dust opacity per dust mass (κ_1.3 mm) of 0.86 cm² g⁻¹ (thick ice mantles coagulated at 10⁵ cm⁻³ from Ossenkopf & Henning 1994) and a gas-to-dust ratio of 100, see Jørgensen et al. (2007). The flux from the unresolved emission is estimated by fitting a point source to baselines longer than 40 kλ, as in Jørgensen et al. (2007), and therefore it avoids contamination from the dense core itself. The result of fitting a point source gives a flux of 27.0 mJy, see Table 9, which implies a mass of

$$\begin{equation} M_{1.3 \,{\rm mm}}=0.024\,{M_{\odot }} \end{equation} \tag{ 5 }$$

for a temperature of 30 K (as used in Jørgensen et al. 2007). This mass will be used as a first estimate for the circumstellar disk mass (see Jørgensen et al. 2007 for a discussion).

The disk-to-dense core mass ratio, M_disk/M_densecore, is estimated using the disk mass from Equation (5) and the dense core mass. This ratio is low (≈0.1) but comparable to Class 0 objects (Enoch et al. 2011; Jørgensen et al. 2009; Enoch et al. 2009a).

A powerful way to constrain the physical parameters of dense cores and YSOs is by fitting the broadband SED (e.g., Robitaille et al. 2007). In the case of L1451-mm, only detections at 160 and 1200 μm are available, which make the broadband SED fit a not well-constrained problem. Here, the information from the 1.3 mm continuum observations (SMA) is extremely important to help discriminate between different physical models (e.g., Enoch et al. 2009a).

In order to compare the continuum emission model with the interferometric observations, the model is sampled in uv-space to match the observations using the uvmodel task in MIRIAD. The synthetic and observed visibilities are both binned in uv-distance (using the uvamp task in MIRIAD), and then they are added as an extra term to the χ² to minimize,

$$\begin{equation} \chi ^2_{{\rm vis}}=\sum _{i}\left(\frac{V_{i,{\rm obs}} - V_{i,{\rm model}}}{\sigma _{V_{i,{\rm obs}}}} \right)^2, \end{equation} \tag{ 6 }$$

where V_{i, obs} and V_{i, model} are the average observed and synthetic visibilities in the i bin, respectively, while $\sigma _{V_{i,{\rm obs}}}$ is the uncertainty of the observed average visibility. The χ² subject to minimization is

$$\begin{equation} \chi ^2 = (1-Q) \chi ^2_{{\rm SED}} + Q \chi ^2_{{\rm vis}}, \end{equation} \tag{ 7 }$$

where Q is an ad hoc weight used to control how important it is to fit the visibilities compared to the broadband SED. In this case a Q = 0.5 is used, which gives the same weight to the SED and visibilities fit.

Because of the problem's high dimensionality, the χ² minimization is carried out with a genetic algorithm (Johnston et al. 2011), while the model SEDs and visibilities are calculated with a new Monte Carlo radiation transfer code (T. P. Robitaille et al. 2011, in preparation), which is based on the radiative transfer code presented by Whitney et al. (2003b). The new code uses ray tracing for the thermal emission at submillimeter and millimeter wavelengths, providing excellent signal to noise to fit the long-wavelength SED and visibilities.

Using this fitting program, we explore three models with increasing levels of complexity to explain our observations: (1) starless isothermal dense core, (2) dense core with a YSO and disk at the center, and (3) dense core with a central FHSC.

The parameter ranges searched using the genetic algorithm are given in Table 10, a summary of the fit results is shown in Figure 9, and the model parameters are listed in Table 11.

Zoom In Zoom Out Reset image size

Table 10. Parameter Ranges Searched

Parameter	Description	Value/Range	Sampling
Starless Dense Core
p	Exponent of density profile	−2–0	Linear
Mdensecore	Envelope mass (M☉)	10−5–10	Logarithmic
Tdust	Envelope temperature (K)	5–30	Linear
Rmax	Envelope outer radius (AU)	1000–5000	Logarithmic
Dense Core with Disk and YSO
Lint	Intrinsic luminosity (L☉)	10−4–0.1	Logarithmic
$\dot{M}_{{\rm infall}}$	Infall rate (M☉ yr−1)	10−9–10−6	Logarithmic
Renv	Outer envelope radius (AU)	1000–5000	Logarithmic
Rcent	centrifugal radius (AU)	50–1000	Logarithmic
Rdisk	outer disk radius (AU)	50–1000	Logarithmic
Mdisk	Disk mass (M☉)	10−8–10−3	Logarithmic
i	Viewing angle (deg)	0–90	Linear

Download table as:  ASCIITypeset image

Table 11. Best Model Fit Parameters

Model	Lint	$\dot{M}_{{\rm infall}}$a	Mdisk	Rdisk	Viewing Angle	Mdensecore	Renv	Tdust
	(L☉)	(M☉ yr−1)	(M☉)	(AU)	(deg)	(M☉)	(AU)	(K)
Starless dense core	...	...	...	...	...	0.353	1001	10.2
Dense core with disk and YSO	0.0450	7.1 × 10−6	0.086	107	64.9	0.146	1007	...
Dense cre with central FHSC	...	...	...	...	...	0.56	2543	...

Note. ^aInfall rate derived using Equation (14) for M_* = 0.086 M_☉.

Download table as:  ASCIITypeset image

Starless isothermal dense core. The simplest model consists of a pure isothermal dense core with a density profile n∝r^−p, where the exponent, p, of the density profile is a free parameter, but constrained to values smaller than 2. The dense core temperature, T_dust, and the outer radius, R_env, are also parameters in the fit, which are listed in Table 11 along with the dense core mass, M_env. The best model fit for the starless isothermal dense core model is shown in red in Figure 9, where it shows that this model does not provide a good match to the visibilities. If the power-law density exponent is not constrained, then a very steep density profile, n∝r^−2.8, can actually match both the SED and visibilities. However, a power-law exponent of −2.8 is beyond the range deemed physically reasonable. Because, even though the outer region of starless cores (and cylindrical filaments) can have similar steep density profiles (e.g., Tafalla et al. 2002), in the inner region (r < 3000 AU) their density profiles are flat, which is exactly the region we are interested in to produce the compact emission.

Dense core with a central YSO and disk. The next model fitted is one composed of a dense core with a central YSO and disk. The dense core is modeled as a rotating and infalling envelope (Ulrich 1976), with outer radius R_env, total mass M_env, and infall rate $\dot{M}_{\rm infall}$. The density of the envelope is given in spherical polar coordinates by

$$\begin{eqnarray} \rho (r,\theta,\phi) &=& \rho _0^{\rm env}\left(\frac{r}{R_{\rm c}}\right)^{-3/2}\left(1 + \frac{\mu }{\mu _0}\right)^{-1/2}\nonumber\\ &&\,\times\,\left(\frac{\mu }{\mu _0} + \frac{2\mu _0^2R_{\rm c}}{r}\right)^{-1}, \end{eqnarray} \tag{ 8 }$$

where R_c is the centrifugal radius, μ = cos θ, and μ₀ is the cosine of the polar angle of a streamline of infalling particles as r → ∞, which is given by

$$\begin{equation} \mu _0^3 + \mu _0\left(\frac{r}{R_{\rm c}} - 1\right) - \mu \left(\frac{r}{R_{\rm c}}\right) = 0. \end{equation} \tag{ 9 }$$

The normalization constant ρ^env₀ is related to the infall rate by

$$\begin{equation} \rho _0^{\rm env} = \frac{\dot{M}_{\rm infall}}{4\pi \left(GM_*R_{\rm c}^3\right)^{1/2}}, \end{equation} \tag{ 10 }$$

where M_* is the mass of the central object. The disk is modeled as a passive flared disk described in cylindrical polar coordinates by

$$\begin{equation} \rho (R,z,\phi) = \rho _0^{\rm disk}\,\left(\frac{100\,{\rm AU}}{R}\right)^{\beta - q}\,\exp {\left[-\frac{1}{2}\left(\frac{z}{h(R)}\right)^2\right]}, \end{equation} \tag{ 11 }$$

where ρ^disk₀ is defined by the disk mass M_disk, q is the surface density radial exponent (which we set to −1), β is the disk flaring power (set to 1.25), and the disk scale height h(R) is given by

$$\begin{equation} h(R) = h_0\,\left(\frac{R}{100\,{\rm AU}}\right)^\beta, \end{equation} \tag{ 12 }$$

where h₀ is the scale height at 100 AU and it is set to 10 AU.

The viewing angle i is a parameter in the fitting. The central protostar is modeled as an object with an effective surface temperature of 3000 K and intrinsic luminosity L_int, where the parameter L_int includes the luminosity due to accretion.

The temperature is computed self-consistently with the density using the radiation transfer code, see Whitney et al. (2003a, 2003b, 2004). It assumes a geometry (e.g., Ulrich envelope model with a flared disk), the dust properties, and local thermodynamic equilibrium. Here we use dust opacities of Ossenkopf & Henning (1994) for dust grains with thick ice mantles after 10⁵ years of coagulation at a density of 10⁶ cm⁻³.

The best model parameters are listed in Table 11. This model provides an excellent fit to the visibilities, while the SED fit underestimates the flux at 160 μm.

Given the best-fit result we estimate the accretion luminosity, L_acc, as

$$\begin{equation} L_{{\rm acc}} = \frac{{\it GM}_{*}\,\dot{M}_{{\rm acc}}}{R_{*}}, \end{equation} \tag{ 13 }$$

where M_* and R_* are the protostellar mass and radius, and $\dot{M}_{{\rm acc}}$ is the accretion rate onto the protostar. Using Equation (13) we obtain an accretion luminosity of 6 L_☉, where we have assumed that the accretion rate is the same as the infall rate,

$$\begin{equation} \dot{M}_{{\rm acc}}=\dot{M}_{{\rm infall}}= 2.4\times 10^{-5}\,\left(\frac{M_*}{{M_{\odot }}}\right)^{0.5}{M_{\odot }\,{\rm yr^{-1}}}, \end{equation} \tag{ 14 }$$

the minimum mass of the central object is the mass of the disk, M_* = M_disk = 0.086 M_☉, and that the central object might be a young protostar, R_* = 3 R_☉. This expected accretion luminosity is 100 times higher than what can be kept hidden at the center of the core in the best-fit model, and therefore some of the assumptions must be clearly misrepresenting reality. Another simple estimate that is derived from Equation (13) is the accretion rate onto the central object needed to produce the same luminosity of the best model, obtaining $\dot{M}_{{\rm acc}}=2.45\times 10^{-8}\,{M_{\odot }\,{\rm yr^{-1}}}$ (much lower than the estimated infall rate, 7.1 × 10⁻⁶ M_☉ yr⁻¹). Clearly, in order to make a self-consistent model an extremely low accretion rate must be assumed.

Dense core with a central FHSC. Another possibility to reconcile the low luminosity observed and the accretion rates expected from the models presented above is to increase the stellar radius, R_*. There is one object which has been predicted from numerical simulations (Larson 1969; Masunaga et al. 1998; Masunaga & Inutsuka 2000; Machida et al. 2008) that is large enough to fit the description: the first hydrostatic core (FHSC).

Masunaga et al. (1998) performed radiation hydrodynamic simulations of a collapsing core until the formation of an FHSC and found that for an initial core of 0.3 M_☉ (model M3a), an FHSC of R_* ∼ 5 AU and L_int ≈ 0.03 L_☉ (including accretion luminosity) is formed. We use the density and temperature profile resulting from this simulation as the input parameter for the radiative transfer code. The mass and radius of the dense core are varied (see Table 11) to find the best match of the observed SED and visibilities, see green points in Figure 9. Note that this model does not use the same density and temperature profile as used in the dense core with disk and YSO model, and this explains the different SEDs. It is clear that the FHSC model visibilities do not provide as good a match to the observations as the YSO model, however, the FHSC model only have two free parameters compared to the seven free parameters of the YSO model.

If there is no source of heating within the core, then ¹²CO should freeze out onto dust grains (e.g., Tafalla et al. 2002). Therefore, the presence of the ¹²CO (2–1) emission in itself is strongly suggestive of a central heating source. From the NH₃, N₂H⁺, and NH₂D centroid velocity maps it is clear that the velocity gradient is in the right ascension (R.A.) direction while the ¹²CO emission is more or less perpendicular. Since the orientation of the dust continuum emission detected with SMA is almost perpendicular to the ¹²CO (2–1) emission, we argue in favor of a central object and outflow system.

The amount of mass needed to keep material at 560 AU with a velocity of 1.3 km s⁻¹ (similar parameters to the ¹²CO emission) is ≈0.53 M_☉, which is almost twice the total mass in the dense core. Therefore, despite the low velocity seen in the ¹²CO (2-1) emission this gas is unbound and it is consistent with ¹²CO (2–1) tracing a slow molecular outflow.

The physical parameters of molecular outflows are typically calculated using both ¹²CO and ¹³CO lines (e.g., Arce et al. 2010). Unfortunately, in L1451-mm there are no detections of ¹³CO (2–1) and we are left to use only the ¹²CO (2–1) emission to study the high-velocity gas.

A lower limit for the mass entrained by the outflow, M_flow, is estimated assuming that the ¹²CO (2–1) emission is optically thin (see details in the Appendix) and with an excitation temperature of T_ex = 20 K. The momentum (P_flow) and energy (E_flow) of the outflow along the line of sight are estimated following Cabrit & Bertout (1990),

$$\begin{eqnarray} P_{{\rm flow}} &=& \sum _i M_{{\rm out},i}\, |v_i -v_{{\rm center}}|, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} E_{{\rm flow}} &=& \frac{1}{2} \sum _i M_{{\rm out},i} \, |v_i -v_{{\rm center}}|^2, \end{eqnarray} \tag{ 16 }$$

where M_{out, i} is the mass in voxel i, v_center is the velocity of the core, and v_i is the velocity of voxel i. Also, the outflow characteristic velocity, v_flow, is calculated as P_flow/M_flow. An upper limit for the lobe size, R_lobe, is estimated from the (not deconvolved) extension of the redshifted and blueshifted ¹²CO emission (2 R_lobe) ≈45, or ≈1120 AU at the distance of Perseus. We also calculate the dynamical time, τ_dyn = R_lobe/v_flow, mechanical luminosity, L_flow = E_flow/τ_dyn, force, F_flow = P_flow/τ_dyn, and rate, $\dot{M}_{{\rm flow}}=M_{{\rm flow}}/\tau _{{\rm dyn}}$.

The outflow properties are calculated using only the voxels with signal-to-noise ratio higher than 2. To avoid contamination in the outflow parameters from the cloud emission the central channel is not used in the calculations, and a second estimate is calculated where the three central channels are removed. The outflow parameters are reported in Table 12, where we note that the differences between the quantities calculated using both methods are small.

The outflow properties presented in Table 12 show an extremely weak outflow in L1451-mm. However, when compared to the outflow found toward L1014 by Bourke et al. (2005), both P_flow and E_flow are close to the lower limits estimated using similar assumptions.

The dynamical time is consistent with the FHSC estimated lifetime (e.g., Machida et al. 2008).

Recent three-dimensional radiation magnetohydrodynamics simulations of the dense core collapse (Machida et al. 2008; Tomida et al. 2010b; Commerçon et al. 2010) show that when the FHSC is formed a slow outflow can be driven even before the existence of a protostar (see also Tomisaka 2002; Banerjee & Pudritz 2006). In these simulations, the outflow that is generated is poorly collimated and typically has maximum velocity of 3 km s⁻¹—very similar to the observed outflow in L1451-mm. In contrast, theoretical models (Shang et al. 2007; Pudritz et al. 2007) and observations (e.g., Arce et al. 2007) indicate that outflows from young protostars are highly collimated and exhibit velocities of a few tens of km s⁻¹ although outflows from VeLLOs display lower velocities (André et al. 1999; Bourke et al. 2005).

The properties of the L1451-mm molecular outflow are consistent with a picture where an FHSC is the driving source of the poorly collimated and slow outflow observed.

If the levels of the molecule are populated following a Boltzmann distribution of temperature T_ex, then the column density can be expressed as

$$\begin{eqnarray} &&N_{J} = \frac{8\pi \nu ^{3}}{c^{3}}\frac{g_{J}}{g_{J+1}} \frac{1}{A_{J+1\rightarrow J}} \frac{\int \tau \, dv}{\left(1-e^{-h\nu /kT_{{\rm ex}}}\right)} N_{J} = 93.28 \left(\frac{2J+1}{2J+3}\right) \frac{\nu ^{3}}{A_{J+1\rightarrow J}}\frac{\int \tau \, dv}{\left(1-e^{-T_{0}/T_{{\rm ex}}}\right)}\,, \end{eqnarray} \tag{ A1 }$$

where g_J = (2J + 1) is the statistical weight of level J for a linear rotor molecule, ν is the transition frequency in units of GHz, A_{J + 1 → J} is the spontaneous emission coefficient in s⁻¹, τ is the transition optical depth, the velocity is in km s⁻¹, and

$$\begin{equation} T_{0}\equiv \frac{h\,\nu }{k}\,. \end{equation} \tag{ A2 }$$

In the case of ¹²CO (2–1), we use ν = 230.538 GHz and A_{2 → 1} = 6.91 × 10⁻⁷ s⁻¹ (obtained from Leiden Atomic and Molecular Database¹¹), and therefore T₀ = 11.0641 K.

Using the equation of radiative transfer to relate the observed emission, T_R, with excitation temperature, T_ex, and background temperature, T_cmb, we obtain

$$\begin{eqnarray} &&T_{R} = T_{0} \left[ \frac{1}{(e^{T_{0}/T_{{\rm ex}}}-1)} - \frac{1}{(e^{T_{0}/T_{{\rm cmb}}}-1)} \right]\left(1-e^{-\tau }\right) \int T_{R}\, dv = T_{0} \left[ \left(e^{T_{0}/T_{{\rm ex}}}-1\right)^{-1} - \left(e^{T_{0}/T_{{\rm cmb}}}-1\right)^{-1} \right] \int \tau \, dv, \end{eqnarray} \tag{ A3 }$$

where optically thin emission is assumed.

Combining Equations (A1) and (A3), the column density of the level J can be calculated as

$$\begin{eqnarray} N_{J} &=& 93.28 \left(\frac{2J+1}{2J+3}\right) \frac{\nu ^{3}}{A_{J+1\rightarrow J}} \frac{\int T_{{\rm MB}}\, dv}{(1-e^{-T_{0}/T_{{\rm ex}}})T_{0} [ (e^{T_{0}/T_{{\rm ex}}}-1)^{-1} - (e^{T_{0}/T_{{\rm cmb}}}-1)^{-1} ]}\nonumber\\ N_{1}({{\rm ^{12}CO}}) &=& 9.924\times 10^{14} \frac{\int T_{{\rm MB}}\, dv}{(1{-}e^{-11.06/T_{{\rm ex}}})11.06 [ (e^{11.06/T_{{\rm ex}}}-1)^{-1} {-} (e^{11.06/T_{{\rm cmb}}}{-}1)^{-1} ]} \ {\rm cm^{-2}}, \end{eqnarray} \tag{ A4 }$$

where Equation (A4) gives the column density of the level J = 1 of ¹²CO using the ¹²CO (2–1) transition emission.

The total column density of ¹²CO is calculated as

$$\begin{eqnarray} &&N({{\rm ^{12}CO}}) = \frac{Z}{g_{J} e^{-E_{J}/k\,T_{{\rm ex}}}} N_{J}({{\rm ^{12}CO}}) = \frac{Z}{(2J+1)}e^{J(J+1) B\, h/k\,T_{{\rm ex}}} N_{J}({{\rm ^{12}CO}}), \end{eqnarray} \tag{ A5 }$$

where B is the rotation constant for a linear rotor (B = 57.635968 GHz for ¹²CO), and Z is the partition function, which can be approximated as Z ≈ k T_ex/(h B). Therefore, in the case of ¹²CO J = 1 we obtain

$$\begin{eqnarray} &&N({{\rm ^{12}CO}}) = 1.2\times 10^{14} \, \frac{T_{{\rm ex}}/11.06 }{(e^{11.06/T_{{\rm ex}}}-1)^{-1} - (e^{11.06/T_{{\rm cmb}}}-1)^{-1}} \frac{e^{5.53/T_{{\rm ex}}} \,\int T_{{\rm MB}}\, dv}{(1-e^{-11.06/T_{{\rm ex}}})} \ {\rm cm^{-2}}, \end{eqnarray} \tag{ A6 }$$

which when combined with a ¹²CO abundance with respect to H₂, [¹²CO/H₂] = 10⁻⁴, provides an estimate of the total column density of H₂.

The final conversion between column density and mass is done using

$$\begin{eqnarray} &&M = 2.71 \times 10^{-8} \left(\frac{N({{\rm ^{12}CO}})}{10^{14}\,{\rm cm^{-2}}}\right) \left(\frac{[{\rm H_2}/{{\rm ^{12}CO}} ]}{10^{4}}\right) \left(\frac{d}{250\,{\rm pc}}\right)^{2} \left(\frac{A_{{\rm sky}}}{\rm arcsec^{2}}\right) \left(\frac{\mu }{2.3}\right) {M_{\odot }} \,, \end{eqnarray} \tag{ A7 }$$

where A_sky is the area on the sky used to calculate N(¹²CO), and μ is the mean molecular weight.