THE SPITZER C2D SURVEY OF NEARBY DENSE CORES. XI. INFRARED AND SUBMILLIMETER OBSERVATIONS OF CB130

A three-color image of CB130-1 using cores2deeper IRAC data is presented in Figure 2, with 3.6 μm as blue, 4.5 μm as green, and 8.0 μm as red. Two YSOs are identified in these data based on their positions in various color–color and color–magnitude diagrams (Harvey et al. 2007). These two YSOs are 15'' apart, which corresponds to 4100 AU at a distance of 270 pc (Straižys et al. 2003). We name the western source (α = 18^h 16^m 164, δ = −02°32'380 [J2000]) CB130-1-IRS1 and the eastern source (α = 18^h 16^m 174, δ = −02°32'411 [J2000]) CB130-1-IRS2 and label these sources in Figure 2.

Zoom In Zoom Out Reset image size

Figure 3 shows a three-color image of CB130-1 using the NIR data with J as blue, H as green, and Ks as red. Both YSOs are detected and labeled in the figure. Cone-shaped nebulosity extends from CB130-1-IRS1 toward CB130-1-IRS2, as marked by arrows in the figure. Similar nebulosity is also observed in L1014 (Huard et al. 2006) and is indirect evidence of an outflow. However, we mapped the CB130-1 region with CO (J = 2 → 1) at the CSO and found no significant evidence of outflowing gas.

Zoom In Zoom Out Reset image size

Figure 4 shows 850 μm contours overlaid on the 1.25–24 μm NIR and Spitzer images. The 850 μm data probes the envelope since dust emission from the envelope becomes optically thin at long wavelengths. The 850 μm contours peak at CB130-1-IRS1 and avoid CB130-1-IRS2. CB130-I-IRS2 is brighter at all wavelengths up to 8 μm, but CB130-1-IRS1 is stronger at 24 μm. These facts support the conclusion by Launhardt et al. (2010) that CB130-1-IRS1 is younger and more embedded than CB130-1-IRS2.

Zoom In Zoom Out Reset image size

A color–color diagram constructed from the IRAC photometry for all sources detected in the Spitzer images of the CB130 region is plotted in Figure 5. Sources inside the box are generally Class II sources (Allen et al. 2004). The two sources in CB130-1 are clearly redder than background stars. The color–color diagram shows that CB130-1-IRS1 and CB130-1-IRS2 are in different phases of evolution. CB130-1-IRS2 is clearly in the region of Class II sources. The sources for which either the [3.6]–[4.5] color is greater than 0.8 mag or the [5.8]–[8.0] color is greater than 1.1 mag are all faint and classified as galaxies except for CB130-1-IRS1 (Harvey et al. 2007). Since no YSO candidates in CB130-2 and CB130-3 were detected with CTIO and Spitzer, we will concentrate on CB130-1 for the remainder of the paper.

Zoom In Zoom Out Reset image size

The photometric data for CB130-1-IRS1 are given in Table 1, which lists wavelength, flux density, uncertainty in flux density, aperture diameter, and telescope. At 3.6–24 μm, we list the flux averaged between the c2d and cores2deeper images, weighted by the inverse of the uncertainties. The spectral energy distribution (SED) is plotted in the upper panel of Figure 6, including the IRS spectrum, which clearly contains 10 μm silicate and 15.2 μm CO₂ absorption features. Using the 15.2 μm CO₂ absorption feature, we can derive the CO₂ abundance toward CB130-1-IRS1. The continuum for the spectrum is constructed by fitting a first-order polynomial to the ranges 12–14.7 and 16.2–18 μm. After subtracting the continuum, we calculated the CO₂ column density from N(CO₂)=∫τ(ν)dν/A, where τ(ν) is an optical depth, ν is the wave number (cm⁻¹), and A = 1.1 × 10⁻¹⁷ (cm) is the band strength (Gerakines et al. 1995). We find ∫τ(ν)dν = 28 ± 0.4 cm⁻¹ and N(CO₂)=(2.6 ± 0.05) × 10¹⁸ cm⁻². We calculated the H₂ column density using 850 μm data. Using Equation (6) in Lee et al. (2003), we have N(H₂)=4.5 × 10²² cm⁻². From N(CO₂) and N(H₂), the abundance of CO₂ ice is 5.8 × 10⁻⁵, which requires a substantial fraction of the carbon.

Zoom In Zoom Out Reset image size

Based on the observed SED, we can calculate standard observational signatures including L_bol, T_bol, and L_bol/L_(>350 μm). We calculate L_bol by integrating the flux over all wavelengths using the average of the results obtained from midpoint and prismoidal integration methods. T_bol is the temperature of a blackbody having the same flux-weighted mean frequency as the source and is calculated following Myers & Ladd (1993). L_(>350 μm) is the luminosity emitted at λ>350 μm. CB130-1-IRS1 has L_bol = 0.23 ± 0.03  L_☉,  T_bol = 59 ± 1.6 K, and L_bol/L_(>350 μm) = 9.1. L_bol/L_(>350 μm) is used to distinguish Class 0 from Class I sources. The pre-protostellar core/Class 0 transition occurs around L_bol/L_(>350 μm) ∼ 35 and the Class 0/I transition occurs at L_bol/L_(>350 μm) ∼ 175 (Young & Evans 2005). Since L_bol/L_(>350 μm) is less than 35, we would guess that CB130-1 is in the pre-protostellar core stage from these criteria, but CB130 clearly contains embedded sources. Another criterion to define a Class 0 object is the bolometric temperature of the source. Class 0 and I objects are divided at T_bol = 70 K (Chen et al. 1995); thus, CB130-1-IRS1 is classified as a Class 0 source based on T_bol.

The bolometric luminosity includes the luminosities of the embedded object, the disk, and the interstellar radiation field (ISRF), which can add several tenths of an L_☉ to L_bol (Evans et al. 2001). The luminosity of the internal source can be estimated in two ways. First, assuming that all emission at λ ⩾ 100 μm arises from external heating, L_int can be calculated by integrating the flux at λ ⩽ 100 μm, giving L_int∼ 0.056 L_☉. Alternatively, we can estimate L_int by applying the relation between the flux at 70 μm, normalized to 140 pc, and the L_int found by Dunham et al. (2008), giving L_int ∼ 0.07 L_☉. These two estimates indicate that CB130-1-IRS1 is in any case a low-luminosity object and may be a VeLLO. However, they are only rough estimates; radiative transfer models of CB130-1 are required to reliably determine the physical properties of the internal source (Section 5).

Assuming that the core is optically thin at submillimeter wavelengths, the envelope mass can be estimated from $M= \frac{S_{\nu } D^2}{B_\nu (T_d) \kappa _\nu }$, where S_ν is the flux density at 850 μm, D is the distance to the object, B_ν(T_d) is the blackbody function, and κ_ν is the opacity of gas and dust at 850 μm (Young et al. 2003). Pontoppidan's dust model (K. M. Pontoppidan et al. 2011, in preparation) gives κ_ν = 1.02 × 10⁻² cm⁻¹ g⁻¹ and κ_ν = 1.8 × 10⁻² cm⁻¹ g⁻¹ for OH5 dust (Ossenkopf & Henning 1994). If we assume the core is isothermal at T_d = 10 K, the estimated envelope mass within a 120'' diameter aperture (radius of 16,200 AU) is 3.5 M_☉ with the Pontoppidan model, and 2.0 M_☉ with the OH5 model. These estimates will be refined later.

The photometric data for CB130-1-IRS2 are given in Table 2. At 3.6–24 μm, we list the flux averaged between the c2d and cores2deeper images, weighted by the inverse of the uncertainties.

The SED of CB130-1-IRS2 is plotted in the lower panel of Figure 6. CB130-1-IRS2 is not detected at λ>24μm. The slope of the SED of CB130-1-IRS2 from 2.17 μm to 24 μm is −0.7, which indicates that CB130-1-IRS2 is a Class II object (Greene et al. 1994). We can estimate the bolometric luminosity of the source by integrating the SED. The bolometric luminosity of CB130-1-IRS2 is L_bol = 0.071 ± 0.003 L_☉ and the bolometric temperature is T_bol = 1150 ± 15 K. These values of L_bol and T_bol are not corrected for extinction along the line of sight and are thus lower limits to the intrinsic values.

Molecular line observations can trace the dynamics in the core, such as infall, rotation, and outflow. The observed molecular lines for CB130-1, CB130-2, and CB130-3 are summarized in Tables 3, 4, and 5, respectively. The detected molecular lines are plotted in Figure 7 as solid lines. Toward CB130-2, the detected molecular lines are plotted in the upper three panels in Figure 8. Toward CB130-3, the detected lines are plotted in the lower three panels in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following section, we will focus on CB130-1-IRS1, since the source is a deeply embedded Class 0 object. We will find the internal luminosity and chemical properties of the source.

We modeled CB130-1-IRS1 with the 1D radiative transfer code, DUSTY (Ivezic et al. 1999). DUSTY calculates the dust temperature profile and SED of a protostar embedded in a dense core. The parameters of the model divide into those associated with the envelope, or core, and those associated with the internal source.

The envelope is characterized by the outer radius, inner radius, density distribution, dust properties, and external radiation field. The outer radius of the envelope is not well constrained by the radiative transfer modeling as long as it is large enough to include all the emission. We fix the outer radius at 35,000 AU, based on the angular size of CB130-1 (30 × 13; Lee & Myers 1999). The inner radius (R_in) will be a free parameter. We assume a power-law radial density distribution for the envelope (n = n₀(r/r₀)^−α), where n₀ is the density at a fiducial radius, r₀. The value of n₀ was constrained to match the long-wavelength emission, and α = 1.8 ±  0.2 was constrained to match the radial intensity profile. See Shirley et al. (2002) for an explanation of this method.

We use the Pontoppidan dust opacity model (K. M. Pontoppidan et al. 2011, in preparation). The envelope is heated by both the embedded source and the ISRF described by Evans et al. (2001). The ISRF can vary in strength and is modified by extinction by a surrounding low-density envelope. To constrain the ISRF, we used an initial model of the envelope and modeled the CO (J = 2 → 1) line, since the CO (J = 2 → 1) line is sensitive to the UV part of the ISRF, and thus the extinction. The best-fit value for the UV part of ISRF is G₀ = 4.5 × 10⁻³, corresponding to A_V = 3. We will discuss this in Section 6. With these dust properties, the density normalization (n₀) yields the total mass of the model envelope within 35,000 AU to be 5.3 M_☉. This mass depends only slightly on the best-fitting parameters for luminosity and inner radius, as long as they are close to the best-fit values. Within 16,200 AU (the size used to estimate the observed mass in Section 4.2), the model contains a mass of 2.2 M_☉, slightly less than the isothermal calculation within that radius.

The internal source is characterized by an internal luminosity (L_int) and a stellar temperature (T_s). The luminosity of the disk is included in the model, but the disk luminosity is not easily separable from the stellar luminosity, so we treat the sum as the free parameter: L_int = L_s + L_disk. The emission from the disk is averaged over all inclination angles and added to the input spectrum (Butner et al. 1994). The 24 μm photometry data constrains the disk parameters. The surface density profile of the model disk follows Σ(r) ∝ r^−p with p = 1.5. We choose a disk temperature distribution following T(r) ∝ r^−q. When the disk is a flat passive disk, q = 0.75. For the flat passive disk, all the luminosity comes from the re-radiation from the central star. A flat passive disk produces less 24 μm emission than observed (lower panel of Figure 9). We varied the index of the temperature distribution until it fitted the observed 24 μm data. The best-fit power-law index of the disk temperature distribution is q = 0.4, which corresponds to a flared disk (Chiang & Goldreich 1997). The disk temperature distribution is also less steep than the flat disk case if the disk has an intrinsic luminosity source, such as accretion and back warming from the envelope in addition to radiation from the central star.

Zoom In Zoom Out Reset image size

DUSTY calculates the radiative transfer from the internal source and the external radiation field. Then, the OBSSPHere (Shirley et al. 2002) program, which convolves the model with a beam, is applied to the DUSTY result to convert models into observed flux densities and model radial profiles at wavelength longer than 100 μm. The model parameters are constrained by the observed SED and the observed radial intensity profile at 850 μm.

We tested a grid of models to obtain the best-fit free parameters. After the constraints described above, the remaining free parameters are T_s, L_int, and R_in. We calculated the reduced χ² to quantify the quality of the fit. The data used to constrain the model covers 3.5 μm ⩽ λ⩽ 850 μm, with a total of nine points of photometry data. We do not include H-band photometry data, because H-band flux comes mainly from scattered photons in the envelope, and DUSTY modeling does not properly treat scattering at short wavelengths. Also, we focus on the photometric data and check the IRS spectrum only for consistency. We have nine data points and three free parameters, so the degree of freedom of reduced χ² is 6.

The χ² plots are in Figure 10. When we plot Figure 10, we hold fixed the other two parameters at the values with the smallest χ² values. We plot in panel (a) of Figure 10 varying T_s from 1000 K to 7000 K. The value of T_s does not critically affect the modeled SED, since all the stellar photons are reprocessed by dust close to the star. The lowest χ² is found at T_s = 3000 ± 500 K. Panel (b) is χ² obtained by varying L_int = L_s + L_disk. The best fit is found at L_int = 0.16 ± 0.005 L_☉. The L_int is a little higher than the VeLLO criterion, but much lower than the accretion luminosity calculated from the mass accretion rate in Shu's (1977) model. Panel (c) is the χ² plot varying R_in. The inner radius of the envelope controls the optical depth. When the optical depth increases, the flux density decreases at short wavelengths because dust absorbs energy at short wavelengths and re-emits at long wavelengths. The best-fit inner radius is R_in = 350 ± 25 AU, which gives τ = 0.22 at 100 μm.

Zoom In Zoom Out Reset image size

The resulting model SED is plotted in Figure 9. The flux density at wavelengths longer than 100 μm comes mainly from the envelope. Since the envelope is bigger than the aperture size, the modeled flux density is higher than the observations. To compare to observations, we used OBSSPHere to convolve the model with the aperture that was used for the photometry. The convolved model data are marked with open circles in Figure 9. The SED shows a good fit to the observed data at 24 μm, 70 μm, 160 μm, 350 μm, and 850 μm. However, the modeled SED shows a 24–30 μm turnover and shallow 15.2 μm CO₂ feature, which is not an ideal fit to the IRS spectrum.

The best-fit result with the 1D radiative transfer model fails to fit the IRS spectrum at wavelength 24–30 μm. The flux density in the mid-infrared (MIR) region mainly comes from the disk emission. There are limitations to models of disks averaged over all inclinations. For a more realistic explanation of the MIR spectrum, we used the 2D axisymmetric Monte Carlo (MC) code, RADMC (Dullemond & Dominik 2004). The RADMC program solves the radiative transfer equation in 2D circumstellar dust configurations. Once a dust temperature distribution is determined, RAYTRACE is used to create smooth SEDs for various inclination angles. As with the 1D models, we use some constraints to determine some parameters before running a fine grid; this reduces the multi-dimensional fitting problem to a manageable size.

We used the same dust opacity and ISRF models used for the 1D model. For the 2D model envelope, we used a slowly rotating cloud core with a small angular momentum (Terebey et al. 1984, hereafter TSC). The initial state of the TSC model is a singular isothermal sphere, with density proportional to r⁻². As the core evolves, the collapse of the isothermal sphere includes the effects of an initially uniform and slow rotation. The density of the TSC envelope is determined by the initial angular velocity, Ω₀, and the time, t. We set Ω₀ = 10⁻¹⁴ s⁻¹, which is the value used in the TSC model. As t increases, the infall radius gets larger and the density distribution gets flattened. The best-fit density distribution is constrained, as described for the 1D model, to correspond to t = 35, 000 years. We use the inner radius and the outer radius as 350 AU and 35,000 AU, respectively, which were determined by the 1D model and the core size. The mass of the envelope that matches the long-wavelength emission is 5.5 M_☉, only slightly larger than the 1D model estimate in Section 5.1.

For the disk density structure, we use the axisymmetric flared disk (Dullemond et al. 2001; Whitney et al. 2003; Pontoppidan et al. 2007b). The central disk structure is as follows:

$$\begin{equation} \rho (R,z)= \frac{\Sigma (R)}{H_p (R) \sqrt{2 \pi }} \rm exp \left[-\frac{z^2}{2H_p(R)^2} \right], \end{equation} \tag{ 1 }$$

where the surface density Σ is

$$\begin{equation} \Sigma (R) = \Sigma _{{\rm disk}} \left(\frac{R}{R_{{\rm disk}}} \right)^{-p} \end{equation} \tag{ 2 }$$

and the disk scale height is given as

$$\begin{equation} \frac{H_p (R)}{R} = \frac{H_{{\rm disk}}}{R_{{\rm disk}}} \left(\frac{R}{R_{{\rm disk}}} \right)^{\alpha _{{\rm fl}}}. \end{equation} \tag{ 3 }$$

We set the disk outer radius as 350 AU, which is the inner radius of the envelope. The disk inner radius is 0.03 AU, set to the location at which a dust temperature is 1500 K. We use most of the other disk model parameters in Whitney et al. (2003). The flaring index α_fl is set to 0.25, and the surface density power-law index is set to p = 1. The disk mass is set as M_disk = 0.01 M_☉.

We ran a grid of models varying three free parameters to fit the observed SED. As for the 1D models, we calculated χ² for the photometric points to determine the best fit. The free parameters are the disk scale height, H_disk/R_disk at R_disk = 100 AU, the internal luminosity, L_int, and the inclination angle. Again, we plot the χ² versus the free parameter, setting the other two free parameters at particular values. Panel (a) of Figure 11 presents χ² with the disk scale height H_disk/R_disk varying in the range between 0.01 (Whitney et al. 2003) and 0.125 (Pontoppidan et al. 2007b). The best fit occurs for H_disk/R_disk = 0.05 ±  0.005. Panel (b) in Figure 11 is the χ² versus internal luminosity. Since the disk is now heated only from the reprocessing of the central star, the internal luminosity affects the SED mostly in the MIR. If the internal luminosity is low, it cannot heat the disk enough. If L_int is high, the modeled SED becomes stronger than the observed flux density everywhere. The best fit is at L_int = 0.14 ± 0.005 L_☉, which is a little bit lower than that found for the 1D model (0.16 L_☉). Panel (c) is the χ² plot as a function of the inclination angle. The inclination angle affects the MIR SED, since the inclination angle affects the disk emission. There is a broad minimum in χ² around an inclination angle of 70°, but the SED for the model at 70° badly underestimates the observations at 24 μm, which are well fitted by a model with an inclination angle of 50°, which however badly overestimates the observations at 70 μm. The SEDs with these two inclination angles are plotted in Figure 12. Neither model matches all of the IRS spectrum either. In particular, the observed CO₂ feature is deeper than produced by any of the models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit parameters are summarized in Table 6. The internal luminosity from the 1D model is 0.16 L_☉, the internal luminosity from the 2D model is 0.14 L_☉. Both luminosities are slightly higher than the VeLLO criterion (L_int = 0.1 L_☉). CB130-1-IRS1 is a low-luminosity source, but not a VeLLO.

We calculate the radial gas temperature distribution using a gas energetics code (Doty & Neufeld 1997; Young et al. 2004b), which includes energy transfer between gas and dust, heating by cosmic rays, photoelectric heating, and molecular cooling.

The collision rate between gas and dust depends on the cumulative grain cross section. In this simulation, the grain cross section per baryon is 6.09 × 10⁻²² cm² (Young et al. 2004b), based on the OH5 model. This is slightly different from the opacity model we have used for DUSTY modeling, but the line modeling does not critically depend on the opacity model. We take the cosmic ray ionization rate as 3.0 × 10⁻¹⁷ s⁻¹ (Roberts & Herbst 2002; van der Tak & van Dishoeck 2000).

The photoelectric heating is determined by the strength of the UV part of the ISRF and the electron number density. The electron number density is assumed to be 0.001 cm⁻³ (Doty et al. 2002). The CO line shows stronger emission when the envelope is heated by the unattenuated ISRF. So we vary the extinction of the ISRF by the surrounding material until it reproduces the observed CO line. See Evans et al. (2005) for a full explanation of this method. The relevant part of the ISRF for the gas heating is in the FUV, where the ISRF is characterized by G₀. Assuming a plane parallel slab, the ISRF is attenuated following $G_0 = G({\rm ISM}) e^{-\tau _{\rm UV}}$, where G(ISM) is the unattenuated FUV field. The best match to the observed CO line requires G₀/G(ISM) = 4.5 × 10⁻³. The relation between optical depth in the UV band and the V band is τ_UV = 1.8A_V. Thus, we attenuate the ISRF at the outer boundary of the envelope by A_V = −0.56ln(4.5 × 10⁻³) = 3. We used this value for the external attenuation of the ISRF in Section 5.1. Of course, the FUV is further attenuated as it propagates into the cloud, and the model accounts for this.

The cooling rate mainly depends on the abundance and the line width of CO. We set the CO fractional abundance relative to H₂ to be 7.4 × 10⁻⁵ (Bacmann et al. 2002; Jørgensen et al. 2002). We use a Doppler b = 0.8 km s⁻¹, which fits best the CO molecular line.

The calculated gas temperature is plotted in the middle panel of Figure 13, along with the dust temperature. The dust temperature and gas temperature are almost the same in the inner region due to the high gas–dust collision rate because of the high density in the inner region, shown in the upper panel. In the outer part of the cloud, the gas temperature is lower than the dust temperature, because the gas temperature has decoupled from the dust temperature and the photoelectric heating from the attenuated ISRF is too low to heat the gas.

Zoom In Zoom Out Reset image size

We use step function abundance profiles. Molecules tend to freeze out onto the grain surfaces in cold and dense regions. The inner part of the cloud has higher density compared to the outer region, so molecules freeze out faster. If the dust temperature becomes high enough, the molecules will sublimate. The CO sublimation temperature is about 20 K. The highest temperature in the model cloud is less than 20 K in the innermost region of the envelope, 0.003 pc, since CB130-1-IRS1 is a low-luminosity source. The dust temperature is too low to sublimate the frozen-out CO. In the outermost part, the density is too low for significant freeze-out. Thus, a step function is a reasonable profile for the cloud. The three parameters of the abundance profile are the undepleted abundance (X₀), depletion radius (r_D), and depletion fraction (f_D). The best-fitting parameters are summarized in Table 7. The abundance profiles are plotted in the last panel in Figure 13.

The line profiles are modeled with the MC code (Choi et al. 1995) to calculate level populations of molecules. The MC code can handle arbitrary 1D distributions of the systematic velocity, the density, the kinetic temperature, the microturbulence, and the abundance self-consistently.

Infalling envelopes sometimes exhibit a stronger blue wing for optically thick lines with a substantial excitation temperature gradient (Evans 1999; Lee et al. 2007). Even though CB130-1 contains a protostar, there is no significant feature of infall in molecular lines. All lines show a single peak, so we do not include any systematic velocity in the models.

Doppler b parameters are determined by the width of each line. We have used b = 0.8 km s⁻¹ for CO. Since the CO line is optically thick and dominated by the outer part of the envelope, the CO line is difficult to fit with the same Doppler b parameter as others. We used b = 0.2 km s⁻¹ for all the other lines, based on the widths of relatively optically thin lines, such as C¹⁸O and H¹³CO.

After MC modeling, a VT program is used to integrate the emission along the line of sight and to match the velocity and spatial resolutions, as well as the beam efficiency. All lines are assumed to be centered at 7.77 km s⁻¹, determined by averaging the LSR velocity of optically thin lines.

By comparing modeled spectra to observations, we can estimate the abundance. The best-fit abundances are given in Table 7. We set the depletion radius as 0.04 pc. We used the isotope ratio for the low-mass protostellar envelope described in Jørgensen et al. (2004). The ratio of O to ¹⁸O is 540, and the ratio between C and ¹³C is 70. CO destroys N₂H⁺ molecules, so N₂H⁺ has a reverse abundance profile compared to the CO abundance. The N₂H⁺ line has a hyperfine structure. The MC code calculates line modeling for each hyperfine line, and then they are added together.

The best-fit results are shown in Figure 7. The modeled C¹⁸O, CS, and N₂H⁺ (J = 1 → 0) lines show reasonable fits to the observed lines, while modeled CO and HCO⁺ lines show a poor fit with the observed lines. The strengths are in a reasonable range for both HCO⁺ and H¹³CO⁺, but the abundance that fits the observed H¹³CO⁺ line produces a self-reversed feature in the modeled HCO⁺ line. The modeled HCO⁺ has a low excitation temperature in the outer region, so it shows absorption features at the center of the line. However, the observed HCO⁺ line does not show a self-reversed structure. The modeled H₂CO also shows a self-reversed structure. The prediction of self-reversed lines that are not observed is a common problem (Carr et al. 1995; Chen et al. 2009), and it suggests that macroturbulent and clumpy three-dimensional models would be more appropriate. Though modeled DCO⁺, CS, and N₂H⁺ (J = 3 → 2) line strengths are weaker than the observed lines, the difference is less than a factor of 2.5.

CB130-1 is a low-luminosity source. There is no significant evidence of a CO outflow, but there is a reflection nebula suggestive of an outflow cavity. CB130-1 could be in the first hydrostatic core (FHSC) stage or in a quiescent stage between episodic accretion bursts. With chemical evolution modeling, we may be able to distinguish these possibilities. Because it takes time to relax the chemistry in the core, evidence of a past high-luminosity stage might be present in the abundances. We use the evolutionary chemo-dynamical model by Lee et al. (2004) to test this idea. This model calculates the chemical evolution of a core from the prestellar core to the embedded protostellar core stage. At each time step, the density profile, the dust temperature, the gas temperature, and abundances are calculated self-consistently.

At the prestellar stage, the density distribution is described by a sequence of Bonnor–Ebert spheres. We adopt the time duration at the preprotostellar core stage in Table 4 in Chen et al. (2009), which has a longer total timescale at the prestellar stage. Once a singular isothermal sphere has been reached with n(r) ∝ r⁻², the inside-out collapse is initiated, and n(r) approaches r^−1.5 inside the infall radius. The infall radius propagates outward at the sound speed through the envelope.

We use three kinds of luminosity evolution models in the accretion stage. First, in Model 1, we adopt the accretion luminosity from Young & Evans (2005) including the FHSC stage (Larson 1969; Boss 1993; Masunaga et al. 1998). The FHSC represents the core in an initial quasi-static contraction stage when it is still molecular. When the core temperature reaches 2000 K, the molecular hydrogen dissociates, so the core has a second collapse and the FHSC stage ends. Models indicate that the FHSC has a short, but poorly determined, lifetime, about 10⁴ years (Boss & Yorke 1995; Masunaga et al. 1998), a very low luminosity, and a high optical depth. While two candidates have recently been proposed (Chen et al. 2010; Enoch et al. 2010), the FHSC has not yet been clearly detected. In Model 1, the FHSC lasts 20,000 years and the collapse of the FHSC occurs in 100 years. When t = 20, 000 years, L_acc ∼ 0.011 L_☉, which is only 10% of the internal luminosity of CB130-1-IRS1. The luminosity increases rapidly after 20,000 years. We adopt a time at the very end of the FHSC stage of Model 1, which has a luminosity much lower than that of CB130-1-IRS1. To match the luminosity of the object, we would have to be catching the source in an extremely short-lived phase while the FHSC is collapsing.

The luminosity evolution of Model 2 is taken from Lee et al. (2004). In Model 2, the FHSC lasts 20,000 years and the transition occurs in 32,000 years, providing a better chance to see a source with the observed luminosity. However, this longer transition stage has no particular theoretical justification. Here, we choose 45,000 years as our desired time since L_acc is 0.16 L_☉ at that time.

The last model of luminosity evolution, Model 3, uses the episodic accretion luminosity model from Dunham et al. (2010a). The episodic accretion is included in a simple, idealized way rather than in a self-consistent model. The model assumes no accretion from disk to star most of the time, but continued infall from envelope to disk. A mass accretion burst occurs when the disk mass reaches 0.2 times the stellar mass. Then mass accretion from the disk to protostar increases to $\dot{M} = 1 \times 10^{-4}$ M_☉yr⁻¹, which is about 100 times the average mass accretion rate. The accretion luminosity depends on $\dot{M}$, so the accretion luminosity also evolves episodically. The envelope mass of CB130-1 is 5.3–5.5 M_☉. So among the three masses in Dunham et al. (2010a), a 3 M_☉ envelope model is the closest. However, the luminosity evolution of the model with a 3 M_☉ envelope never falls below 1 L_☉ after the FHSC stage ends, so we take the luminosity evolution of the 1 M_☉ envelope model in Dunham et al. (2010a). In Model 3, we adopt 78,000 years as our desired time step. At 78,000 years, the internal luminosity is 0.14 L_☉. The cloud has gone through six burst episodes, and 3000 years have passed after the sixth burst. The luminosity evolution of the three models are plotted in Figure 14. The early stage evolution is the same in all three models. After 20,000 years, Model 2 has the lowest luminosity among the three models.

Zoom In Zoom Out Reset image size

We calculate the dust temperature distribution with DUSTY at each time step. We use the same parameters described in Section 5.1. Then we calculate the gas temperature considering the parameters described in Section 6.1. The chemical evolution is calculated for each of the 512 gas parcels as it falls into the central region. Gas inside the infall radius carries the memory of the conditions from farther out. The chemical calculation includes interactions between gas and dust grains and gas-phase reactions, but does not include chemical reactions on grain mantles explicitly. We assume the surface binding energy of species onto bare SiO₂ dust grains. We use an updated chemical network with a new binding energy of N₂ and deuterated species, which is used in Chen et al. (2009).

The abundance profile at each time step is calculated through the chemical evolution model. We use the same isotope ratio as in Section 6.3. The abundance profiles are plotted in Figure 15. The abundance profiles of Model 1 and Model 2 do not have a CO evaporation region at the center. In both cases, the central source has a low luminosity, and it does not heat the envelope enough to make CO sublimate. Model 2 has higher abundances compared to Model 1, except for DCO⁺. In Model 3, the abundances of CO and C¹⁸O show the CO sublimation region at the central region because of past periods of high luminosity. The modeled abundances of HCO⁺, H¹³CO⁺ are all low compared to the best-fit step function abundances.

Zoom In Zoom Out Reset image size

Once all the physical parameters are determined, we use MC code (Choi et al. 1995) to calculate the level populations. Then we use VT code to produce line profiles. The molecular line profiles with Model 1 are plotted in Figure 16, and with Model 2 in Figure 17. The modeled lines of HCO⁺, H¹³CO⁺, CS, and N₂H⁺ are weaker than the observed lines, and the DCO⁺ line is stronger than the observed line. N₂H⁺ shows a better fit with Model 2, but the J = 3 → 2 line is still weak. The molecular lines with Model 3 are plotted in Figure 18. In Model 3, the C¹⁸O line is much stronger than observed lines. On the other hand, HCO⁺, H¹³CO⁺, and N₂H⁺ are weaker than observed lines. The modeled N₂H⁺ (3–2) line is very weak. All three models show strong C¹⁸O and weak N₂H⁺ lines compared to the observed lines. The modeled C¹⁸O is much stronger in Model 3, because the episodic bursts have sublimated the CO in the inner region. The N₂H⁺ line is weak in all three models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All the models have caveats in their luminosity evolution scenarios. Model 1 does not actually match the observed L_int ∼ 0.14–0.16 L_☉. Model 2 assumes that the transition stage of FHSC lasts for 32,000 years, which is unrealistic. It can give a source luminosity of 0.16 L_☉, and gives a slightly better fit to the C¹⁸O and N₂H⁺ compared to Model 1. However, the modeled DCO⁺ line is too strong, and the modeled CS line is too weak. In Model 3, the envelope mass is lower than our estimates for CB130-1. Also, the best-fit TSC model density from the 2D model is 35,000 years, while the desired luminosity of Model 3 is reached at t = 78, 000 years, which means that more envelope material from Model 3 has already fallen in. The problems with the luminosity evolution scenarios show how chemical modeling can provide additional constraints that are distinct from those from modeling the continuum.

We have one more trick to try. As noted in Section 4.2, a substantial amount of carbon is contained in CO₂ ice. Once CO freezes out onto dust grain surfaces, a fraction of CO ice turns into CO₂ ice (Pontoppidan et al. 2007a). CO₂ ice does not sublimate during the accretion burst even when the CO evaporation temperature is reached. So the amount of CO returned to the gas phase is less than the amount of CO frozen onto grain surfaces during the quiescent phase. The net effect is to decrease gas phase CO and enhance N₂H⁺. We used the same episodic accretion luminosity evolution model as Model 3, and varied the fraction of CO ice which turns into CO₂ ice in the chemical network (Model 4). When 80% of CO ice turns into CO₂ ice, the modeled C¹⁸O and N₂H⁺ lines match well with the observed lines. Also, DCO⁺ and CS fit the data better than in Model 3. The abundance of this model is plotted in Figure 15 and the molecular line profiles are plotted in Figure 19. The CO₂ ice column density obtained from the best-fit model is 2.94 × 10¹⁸ cm⁻², 1.1 times observed value, within the likely range of systematic uncertainties.

Zoom In Zoom Out Reset image size

Among the models, the episodic model including CO₂ ice formation fits best to the observed lines in all molecules. Based on the modeling result, we suggest that CB130-1 has gone through episodic accretion, forming CO₂ ice out of CO ice in quiescent phases. This model can explain the low luminosity of the source, the strong CO₂ ice feature, and most of the gas phase emission lines. There are still inconsistencies in this model, so further development of episodic accretion models is needed.